Технологический практикум: самые продаваемые книги¶

Гореленкова Анастасия, группа 316¶

Подготовка к работе¶

Разрешаем совместную работу на Python и R:

In [1]:
%load_ext rpy2.ipython

Устанавливаем необходимые библиотеки Python:

In [2]:
import pip

import warnings
warnings.simplefilter('ignore')
warnings.filterwarnings('ignore') # предупреждающие сообщения не выводятся



pip.main(['install', 'pandas'])
import pandas as pd

pip.main(['install', 'matplotlib'])
import matplotlib
import matplotlib.pyplot as plt
#графики отображаются сразу после кода
%matplotlib inline

pip.main(['install', 'numpy'])
import numpy as np

pip.main(['install', 'seaborn'])
import seaborn as sns
sns.set(style='darkgrid')

pip.main(['install', 'math'])
import math



pip.main(['install', 'ast'])
import ast # literal_eval

pip.main(['install', 'outlier_utils'])
from outliers import smirnov_grubbs as grubbs # критерий Граббса

pip.main(['install', 'statsmodels'])
import statsmodels.api as sm # qq-график
from statsmodels.stats.diagnostic import lilliefors # критерий Колмогорова-Смирнова в мод. Лиллиефорса
from statsmodels.stats.contingency_tables import mcnemar # тест МакНемара
from statsmodels.stats.outliers_influence import variance_inflation_factor # фактор инфляции дисперсии
from statsmodels.tools.tools import add_constant
from statsmodels.formula.api import ols # многофакторный дисперсионный анализ

pip.main(['install', 'scipy'])
import scipy.stats as stats
from scipy.stats import kstest # критерий Колмогорова-Смирнова
from scipy.stats import shapiro # критерий Шапиро_Уилка
from scipy.stats import anderson # критерий Андерсона-Дарлинга
from scipy.stats import anderson_ksamp
from scipy.stats import cramervonmises # критерий Крамера фон Мизеса
from scipy.stats import cramervonmises_2samp
from scipy.stats import mannwhitneyu # критерий Уилкоксона-Манна-Уитни
from scipy.stats import levene # критерий Левене
from scipy.stats import bartlett # критерий Баретта
from scipy.stats import fligner # критерий Флигнера-Килина
from scipy.stats import pearsonr # коэффициент корреляции Пирсона
from scipy.stats import spearmanr # коэффициент корреляции Спирмена
from scipy.stats import kendalltau # коэффициент корреляции Кендалла
from scipy.stats import chi2_contingency # критерий хи-квадрат
from scipy.stats import fisher_exact # точный тест Фишера
from scipy.stats import f_oneway # однофакторный дисперсионный анализ

pip.main(['install', 'sfrancia'])
from sfrancia import shapiroFrancia # критерий Шапиро-Франсия

pip.main(['install', 'statistics'])
import statistics # mean и т.д.

pip.main(['install', 'pingouin'])
from pingouin import ttest # критерий Стьюдента

pip.main(['install', 'cmh'])
from cmh import CMH # критерий Кохрана-Мантеля-Хензеля
/usr/local/lib/python3.10/dist-packages/_distutils_hack/__init__.py:33: UserWarning: Setuptools is replacing distutils.
  warnings.warn("Setuptools is replacing distutils.")
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: pandas in /usr/local/lib/python3.10/dist-packages (1.5.3)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas) (2023.3.post1)
Requirement already satisfied: numpy>=1.21.0 in /usr/local/lib/python3.10/dist-packages (from pandas) (1.23.5)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.8.1->pandas) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: matplotlib in /usr/local/lib/python3.10/dist-packages (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (4.45.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.4.5)
Requirement already satisfied: numpy>=1.20 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (1.23.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (23.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib) (2.8.2)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (1.23.5)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: seaborn in /usr/local/lib/python3.10/dist-packages (0.12.2)
Requirement already satisfied: numpy!=1.24.0,>=1.17 in /usr/local/lib/python3.10/dist-packages (from seaborn) (1.23.5)
Requirement already satisfied: pandas>=0.25 in /usr/local/lib/python3.10/dist-packages (from seaborn) (1.5.3)
Requirement already satisfied: matplotlib!=3.6.1,>=3.1 in /usr/local/lib/python3.10/dist-packages (from seaborn) (3.7.1)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (4.45.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (23.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib!=3.6.1,>=3.1->seaborn) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=0.25->seaborn) (2023.3.post1)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->matplotlib!=3.6.1,>=3.1->seaborn) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
ERROR: Could not find a version that satisfies the requirement math (from versions: none)
ERROR: No matching distribution found for math
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Collecting ast
  Using cached AST-0.0.2.tar.gz (19 kB)
  Preparing metadata (setup.py): started
  error: subprocess-exited-with-error
  
  × python setup.py egg_info did not run successfully.
  │ exit code: 1
  ╰─> See above for output.
  
  note: This error originates from a subprocess, and is likely not a problem with pip.
  Preparing metadata (setup.py): finished with status 'error'
error: metadata-generation-failed

× Encountered error while generating package metadata.
╰─> See above for output.

note: This is an issue with the package mentioned above, not pip.
hint: See above for details.
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: outlier_utils in /usr/local/lib/python3.10/dist-packages (0.0.5)
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from outlier_utils) (1.23.5)
Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (from outlier_utils) (1.11.4)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: statsmodels in /usr/local/lib/python3.10/dist-packages (0.14.0)
Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (1.23.5)
Requirement already satisfied: scipy!=1.9.2,>=1.4 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (1.11.4)
Requirement already satisfied: pandas>=1.0 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (1.5.3)
Requirement already satisfied: patsy>=0.5.2 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (0.5.3)
Requirement already satisfied: packaging>=21.3 in /usr/local/lib/python3.10/dist-packages (from statsmodels) (23.2)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.0->statsmodels) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.0->statsmodels) (2023.3.post1)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from patsy>=0.5.2->statsmodels) (1.16.0)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (1.11.4)
Requirement already satisfied: numpy<1.28.0,>=1.21.6 in /usr/local/lib/python3.10/dist-packages (from scipy) (1.23.5)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: sfrancia in /usr/local/lib/python3.10/dist-packages (1.0.8)
Requirement already satisfied: numpy in /usr/local/lib/python3.10/dist-packages (from sfrancia) (1.23.5)
Requirement already satisfied: scipy in /usr/local/lib/python3.10/dist-packages (from sfrancia) (1.11.4)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: statistics in /usr/local/lib/python3.10/dist-packages (1.0.3.5)
Requirement already satisfied: docutils>=0.3 in /usr/local/lib/python3.10/dist-packages (from statistics) (0.18.1)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: pingouin in /usr/local/lib/python3.10/dist-packages (0.5.3)
Requirement already satisfied: numpy>=1.19 in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.23.5)
Requirement already satisfied: scipy>=1.7 in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.11.4)
Requirement already satisfied: pandas>=1.0 in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.5.3)
Requirement already satisfied: matplotlib>=3.0.2 in /usr/local/lib/python3.10/dist-packages (from pingouin) (3.7.1)
Requirement already satisfied: seaborn>=0.11 in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.12.2)
Requirement already satisfied: statsmodels>=0.13 in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.14.0)
Requirement already satisfied: scikit-learn in /usr/local/lib/python3.10/dist-packages (from pingouin) (1.2.2)
Requirement already satisfied: pandas-flavor>=0.2.0 in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.6.0)
Requirement already satisfied: outdated in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.2.2)
Requirement already satisfied: tabulate in /usr/local/lib/python3.10/dist-packages (from pingouin) (0.9.0)
Requirement already satisfied: contourpy>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (1.2.0)
Requirement already satisfied: cycler>=0.10 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (0.12.1)
Requirement already satisfied: fonttools>=4.22.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (4.45.1)
Requirement already satisfied: kiwisolver>=1.0.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (1.4.5)
Requirement already satisfied: packaging>=20.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (23.2)
Requirement already satisfied: pillow>=6.2.0 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (9.4.0)
Requirement already satisfied: pyparsing>=2.3.1 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (3.1.1)
Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from matplotlib>=3.0.2->pingouin) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1.0->pingouin) (2023.3.post1)
Requirement already satisfied: xarray in /usr/local/lib/python3.10/dist-packages (from pandas-flavor>=0.2.0->pingouin) (2023.7.0)
Requirement already satisfied: patsy>=0.5.2 in /usr/local/lib/python3.10/dist-packages (from statsmodels>=0.13->pingouin) (0.5.3)
Requirement already satisfied: setuptools>=44 in /usr/local/lib/python3.10/dist-packages (from outdated->pingouin) (67.7.2)
Requirement already satisfied: littleutils in /usr/local/lib/python3.10/dist-packages (from outdated->pingouin) (0.2.2)
Requirement already satisfied: requests in /usr/local/lib/python3.10/dist-packages (from outdated->pingouin) (2.31.0)
Requirement already satisfied: joblib>=1.1.1 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->pingouin) (1.3.2)
Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->pingouin) (3.2.0)
Requirement already satisfied: six in /usr/local/lib/python3.10/dist-packages (from patsy>=0.5.2->statsmodels>=0.13->pingouin) (1.16.0)
Requirement already satisfied: charset-normalizer<4,>=2 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (3.3.2)
Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (3.6)
Requirement already satisfied: urllib3<3,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (2.0.7)
Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests->outdated->pingouin) (2023.11.17)
WARNING: pip is being invoked by an old script wrapper. This will fail in a future version of pip.
Please see https://github.com/pypa/pip/issues/5599 for advice on fixing the underlying issue.
To avoid this problem you can invoke Python with '-m pip' instead of running pip directly.
Requirement already satisfied: cmh in /usr/local/lib/python3.10/dist-packages (1.0.1)
Requirement already satisfied: numpy>=1.18 in /usr/local/lib/python3.10/dist-packages (from cmh) (1.23.5)
Requirement already satisfied: pandas>=1 in /usr/local/lib/python3.10/dist-packages (from cmh) (1.5.3)
Requirement already satisfied: scikit-learn in /usr/local/lib/python3.10/dist-packages (from cmh) (1.2.2)
Requirement already satisfied: python-dateutil>=2.8.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1->cmh) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in /usr/local/lib/python3.10/dist-packages (from pandas>=1->cmh) (2023.3.post1)
Requirement already satisfied: scipy>=1.3.2 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->cmh) (1.11.4)
Requirement already satisfied: joblib>=1.1.1 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->cmh) (1.3.2)
Requirement already satisfied: threadpoolctl>=2.0.0 in /usr/local/lib/python3.10/dist-packages (from scikit-learn->cmh) (3.2.0)
Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.8.1->pandas>=1->cmh) (1.16.0)

Скачиваем функции, которые понадобятся для работы на R:

In [3]:
%%R

oldw <- getOption("warn") # предупреждающие сообщения не выводятся
options(warn = -1)

install.packages("outliers") # поиск выбросов (критерии Граббса и Диксона)
library(outliers)
install.packages("nortest") # критерии Андерсона-Дарлинга, Лиллиефорса, Шапиро-Франсия; Пирсона
library(nortest)
install.packages("goftest") # критерий Крамера фон Мизеса
library(goftest)

install.packages("stats") # rnorm, qq-график
library(stats)
install.packages("boot") # метод огибающих
library(boot)

install.packages("corrplot") # корреляционная матрица
library(corrplot)
install.packages("car") # vif
library(car)

install.packages("data.table") # %like%
library(data.table)
install.packages("dplyr") # %>%
library(dplyr)
install.packages("relevance") # dropNA
library(relevance)
WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/outliers_0.15.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 15967 bytes (15 KB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 15 KB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/nortest_1.0-4.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 6179 bytes

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 6179 bytes


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/goftest_1.2-3.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 32319 bytes (31 KB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 31 KB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: 
Attaching package: ‘goftest’


WARNING: R[write to console]: The following objects are masked from ‘package:nortest’:

    ad.test, cvm.test


WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/boot_1.3-28.1.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 236854 bytes (231 KB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 231 KB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/corrplot_0.92.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 3765850 bytes (3.6 MB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 3.6 MB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: corrplot 0.92 loaded

WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/car_3.1-2.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 579829 bytes (566 KB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 566 KB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: Loading required package: carData

WARNING: R[write to console]: 
Attaching package: ‘car’


WARNING: R[write to console]: The following object is masked from ‘package:boot’:

    logit


WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/data.table_1.14.10.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 5311651 bytes (5.1 MB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 5.1 MB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: data.table 1.14.10 using 1 threads (see ?getDTthreads).  Latest news: r-datatable.com

WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/dplyr_1.1.4.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 1207521 bytes (1.2 MB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 1.2 MB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: 
Attaching package: ‘dplyr’


WARNING: R[write to console]: The following objects are masked from ‘package:data.table’:

    between, first, last


WARNING: R[write to console]: The following object is masked from ‘package:car’:

    recode


WARNING: R[write to console]: The following objects are masked from ‘package:stats’:

    filter, lag


WARNING: R[write to console]: The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


WARNING: R[write to console]: Installing package into ‘/usr/local/lib/R/site-library’
(as ‘lib’ is unspecified)

WARNING: R[write to console]: trying URL 'https://cran.rstudio.com/src/contrib/relevance_2.0.tar.gz'

WARNING: R[write to console]: Content type 'application/x-gzip'
WARNING: R[write to console]:  length 889295 bytes (868 KB)

WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: =
WARNING: R[write to console]: 

WARNING: R[write to console]: downloaded 868 KB


WARNING: R[write to console]: 

WARNING: R[write to console]: 
WARNING: R[write to console]: The downloaded source packages are in
        ‘/tmp/RtmpbiptEF/downloaded_packages’
WARNING: R[write to console]: 
WARNING: R[write to console]: 

WARNING: R[write to console]: 
Attaching package: ‘relevance’


WARNING: R[write to console]: The following object is masked from ‘package:dplyr’:

    last


WARNING: R[write to console]: The following object is masked from ‘package:data.table’:

    last


Особенности данных¶

В работе (в различных заданиях) использованы 2 датасета.

Первый набор данных основан на датасете best-selling-books. В нём содержится информация о самых продаваемых книгах (бестселлерах[1](#01)):

  • название (столбец title)
  • автор(ы) (author)
  • язык оригинала (language)
  • год первой публикации (first_published)
  • приблизительный объём продаж, млн (sales)
  • жанры (genres)

Во-первых, для удобства обработки произведения, написанные в соавторстве, разделены на строки, соответствующие отдельным авторам.

Во-вторых, данные были дополнены следующей информацией, взятой из открытых источников:

  • пол автора (gender)
  • год рождения (birth_year)
  • возраст на момент первой публикации (age), который вычислялся как разница между годом первой публикации публикации и годом рождения

В-третьих, из-за большого количества пропусков и запутанных обозначений был полностью изменён столбец genres. Для единообразия обозначений (в том числе со вторым датасетом) использовались данные сайта Goodreads.


1: Под ними понимаются те книги, которые имеют наибольшее ожидаемое количество проданных экземпляров, а не количество экземпляров, напечатанных или находящихся в собственности в настоящее время; также в список не входят комиксы и учебная литература.

Получившийся набор данных выглядит следующим образом:

Python:

In [4]:
df1 = pd.read_csv('best-selling-books-upgrade.csv')

# заполняем пропуски в genres
df1.genres = df1.genres.fillna('[]')
# преобразуем данные столбца genres из строк в списки
df1.genres = df1.genres.apply(ast.literal_eval)

print('Количество строк:', df1.shape[0])
df1.head(5)
Количество строк: 175
Out[4]:
title author gender birth_year age language first_published sales genres
0 And Then There Were None Agatha Christie F 1890 49 English 1939 100.0 [Classics, Fiction, Crime, Thriller, Mystery T...
1 The Plague (La Peste) Albert Camus M 1913 34 French 1947 12.0 [Fiction, Classics, Philosophy, France, Litera...
2 The Stranger (L'Étranger) Albert Camus M 1913 29 French 1942 10.0 [Classics, Fiction, Philosophy, France, Litera...
3 The Joy of Sex Alex Comfort M 1920 52 English 1972 10.0 [Nonfiction, Sexuality, Relationships, Self He...
4 The Young Guard (Молодая гвардия) Alexander Alexandrovich Fadeyev M 1901 44 Russian 1945 26.0 [Russia, Fiction, Russian Literature, Historic...

R:

In [5]:
%%R

df1 <- read.csv(file = "best-selling-books-upgrade.csv", header = TRUE)

# заполняем пропуски в genres
df1$genres[is.na(df1$genres)] <- "[]"
# преобразуем данные столбца genres из строк в списки
df1$genres <- gsub("\\[|\\]|\\'|\\ ", "", df1$genres)
df1$genres <- strsplit(df1$genres, "\\,")

cat("Количество строк:", nrow(df1), "\n")
head(df1, n=5)
Количество строк: 175 
                              title                          author gender
1          And Then There Were None                 Agatha Christie      F
2             The Plague (La Peste)                    Albert Camus      M
3         The Stranger (L'Étranger)                    Albert Camus      M
4                    The Joy of Sex                    Alex Comfort      M
5 The Young Guard (Молодая гвардия) Alexander Alexandrovich Fadeyev      M
  birth_year age language first_published sales
1       1890  49  English            1939   100
2       1913  34   French            1947    12
3       1913  29   French            1942    10
4       1920  52  English            1972    10
5       1901  44  Russian            1945    26
                                                                                                           genres
1                                         Classics, Fiction, Crime, Thriller, MysteryThriller, Audiobook, Mystery
2                                     Fiction, Classics, Philosophy, France, Literature, FrenchLiterature, Novels
3                                     Classics, Fiction, Philosophy, France, Literature, Novels, FrenchLiterature
4                                  Nonfiction, Sexuality, Relationships, SelfHelp, Reference, Erotica, Psychology
5 Russia, Fiction, RussianLiterature, HistoricalFiction, War, Classics, WorldWarII, Novels, Unfinished, Childrens

Второй использованный датасет - Best Books (10k) Multi-Genre Data. В нём содержится информация о 10,000 самых рекоммендованных книгах согласно рейтингу сайта Goodreads. В работе используются следующие данные:

  • название (title)
  • автор (author)
  • средняя оценка пользователей от 1 до 5 (average_rating)
  • количество оценок (number_of_ratings)
  • жанры (genres)

Так выглядит этот набор данных:

Python:

In [6]:
df2 = pd.read_csv('goodreads-data-upgrade.csv')

# заполняем пропуски в genres
df2.genres = df2.genres.fillna('[]')
# преобразуем данные столбца genres из строк в списки
df2.genres = df2.genres.apply(ast.literal_eval)

# преобразуем данные столбца number_of_ratings из строк в числа
df2['number_of_ratings'] = df2['number_of_ratings'].apply(lambda x: int(x.replace(',', '')))

print('Количество строк:', df2.shape[0])
df2.head(5)
Количество строк: 10000
Out[6]:
title author average_rating number_of_ratings genres
0 To Kill a Mockingbird Harper Lee 4.27 5691311 [Classics, Fiction, Historical Fiction, School...
1 Harry Potter and the Philosopher’s Stone (Harr... J.K. Rowling 4.47 9278135 [Fantasy, Fiction, Young Adult, Magic, Childre...
2 Pride and Prejudice Jane Austen 4.28 3944155 [Classics, Fiction, Romance, Historical Fictio...
3 The Diary of a Young Girl Anne Frank 4.18 3488438 [Classics, Nonfiction, History, Biography, Mem...
4 Animal Farm George Orwell 3.98 3575172 [Classics, Fiction, Dystopia, Fantasy, Politic...

R:

In [7]:
%%R

df2 <- read.csv(file = "goodreads-data-upgrade.csv", header = TRUE)

# заполняем пропуски в genres
df2$genres[is.na(df2$genres)] <- "[]"
# преобразуем данные столбца genres из строк в списки
df2$genres <- gsub("\\[|\\]|\\'|\\ ", "", df2$genres)
df2$genres <- strsplit(df2$genres, "\\,")

# преобразуем данные столбца number_of_ratings из строк в числа
df2$number_of_ratings <- as.integer(gsub(",", "", df2$number_of_ratings))

cat("Количество строк:", nrow(df2), "\n")
head(df2, n=5)
Количество строк: 10000 
                                                        title        author
1                                       To Kill a Mockingbird    Harper Lee
2 Harry Potter and the Philosopher’s Stone (Harry Potter, #1)  J.K. Rowling
3                                         Pride and Prejudice   Jane Austen
4                                   The Diary of a Young Girl    Anne Frank
5                                                 Animal Farm George Orwell
  average_rating number_of_ratings
1           4.27           5691311
2           4.47           9278135
3           4.28           3944155
4           4.18           3488438
5           3.98           3575172
                                                                            genres
1 Classics, Fiction, HistoricalFiction, School, Literature, YoungAdult, Historical
2            Fantasy, Fiction, YoungAdult, Magic, Childrens, MiddleGrade, Classics
3 Classics, Fiction, Romance, HistoricalFiction, Literature, Historical, Audiobook
4          Classics, Nonfiction, History, Biography, Memoir, Historical, Holocaust
5               Classics, Fiction, Dystopia, Fantasy, Politics, School, Literature

1. Апроксимация распределений с помощью ядерных оценок плотности¶

Ядерные оценки позволяют строить гладкие оценки плотностей распределения. Часто их изображают вместе с соответствующими гистограммами. Рассмотрим, например, распределение годов первой публикации бестселлеров.


Для определения количества интервалов разбиения hist() в R по умолчанию использует метод Стёрджеса , а histplot() в Python - максимум оценок Стёрджеса и Фридмана-Диакониса. Метод Стёрджеса в данном случае даёт очевидно неудовлетворительный результат, поэтому в обоих случаях использован метод Фридмана-Диакониса.

Тем не менее, отличия в графиках присутствуют. Видимо, метод по-разному реализован на R и Python: границы отрезков разбиения и их количество не совпадают.

Полоса пропускания ядра установлена по умолчанию.

Python:

In [8]:
edges = np.histogram_bin_edges(df1['first_published'], bins='fd')

print('Количество отрезков разбиения:', edges.shape[0] - 1)
edges
Количество отрезков разбиения: 42
Out[8]:
array([1304., 1321., 1338., 1355., 1372., 1389., 1406., 1423., 1440.,
       1457., 1474., 1491., 1508., 1525., 1542., 1559., 1576., 1593.,
       1610., 1627., 1644., 1661., 1678., 1695., 1712., 1729., 1746.,
       1763., 1780., 1797., 1814., 1831., 1848., 1865., 1882., 1899.,
       1916., 1933., 1950., 1967., 1984., 2001., 2018.])
In [9]:
plt.figure(figsize=(10, 6))

sns.histplot(data=df1['first_published'], stat='density', bins=edges, kde=True)

plt.xticks(fontsize=9)
plt.xlabel('Год', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('Распределение годов первой публикации', fontsize=10, fontweight='bold')
plt.show()

R:

In [10]:
%%R

edges = hist(df1$first_published, freq = FALSE, breaks = "FD", plot = FALSE)

cat("Количество отрезков разбиения:", length(edges$breaks) - 1, "\n")
edges$breaks
Количество отрезков разбиения: 36 
 [1] 1300 1320 1340 1360 1380 1400 1420 1440 1460 1480 1500 1520 1540 1560 1580
[16] 1600 1620 1640 1660 1680 1700 1720 1740 1760 1780 1800 1820 1840 1860 1880
[31] 1900 1920 1940 1960 1980 2000 2020
In [11]:
%%R -w 1000 -h 600

hist(df1$first_published, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Год", ylab = "",
     main = "Распределение годов первой публикации")

lines(density(df1$first_published, adjust=), col = "royalblue", lwd = 2)

Видно, что пик распределения приходится примерно на 80е - 90е годы. Действительно это время, когда, с одной стороны, книги активно пишутся и издаются и уже давно доступны большому количеству людей, но, с другой стороны, ещё не распространён интернет и т.д.

Если строить ядерные оценки плотностей для разных групп, то можно сравнивать распределения в них. Рассмортим, например, как отличаются распределения для авторов разного пола.


Признак gender выбран намеренно, так как выделенные с помощью него подгруппы не пересекаются и любая книга обязательно попадает в одну из подгрупп. В противном случее использовать нормирование (площадь под гистограммой равна 1) было бы неразумно (например, из-за того, что большинство книг имеют не один жанр, некоторые из них при построении отдельных графиков могут учитываться несколько раз и т.д.).

Интервалы разбиения у графиков такие же, как и у графиков выше, полоса пропускания по умолчанию.

Python:

In [12]:
plt.figure(figsize=(10, 6))

sns.histplot(data=df1[df1.gender == 'F'].first_published, stat='density', kde=True,
             bins=edges,
             color='red', alpha=0.3,
             label='Female')
sns.histplot(data=df1[df1.gender == 'M'].first_published, stat='density', kde=True,
             bins=edges,
             color='blue', alpha=0.3,
             label='Male')

plt.xticks(fontsize=9)
plt.xlabel('Год', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('Распределение годов первой публикации', fontsize=10, fontweight='bold')
plt.legend(title='Пол', title_fontsize=10, fontsize=10)
plt.show()

R:

In [13]:
%%R
# функция, которая создаёт прозрачные версии цветов

t_col <- function(color, percent = 50, name = NULL) {
  #      color = color name
  #    percent = % transparency
  #       name = an optional name for the color

# Get RGB values for named color
rgb.val <- col2rgb(color)

# Make new color using input color as base and alpha set by transparency
t.col <- rgb(rgb.val[1], rgb.val[2], rgb.val[3],
             max = 255,
             alpha = (100 - percent) * 255 / 100,
             names = name)

# Save the color
invisible(t.col)
}
In [14]:
%%R -w 1000 -h 600

breaks = hist(df1$first_published, breaks = 145, plot = FALSE)$breaks

my_red <- t_col("red", perc = 60)
hist(df1[df1$gender %like% "F", "first_published"], freq = FALSE,
     breaks = edges$breaks,
     xlim = c(1300, 2020),
     col = my_red,
     xlab = "Год", ylab = "",
     main = "Распределение годов первой публикации")

my_blue <- t_col("royalblue", perc = 60)
hist(df1[df1$gender %like% "M", "first_published"], freq = FALSE,
     breaks = edges$breaks,
     col = my_blue,
     add = TRUE)

lines(density(df1[df1$gender %like% "F", "first_published"], adjust=), col = "red", lwd = 3)
lines(density(df1[df1$gender %like% "M", "first_published"], adjust=), col = "royalblue", lwd = 3)

legend("topleft", inset=c(0.03, 0), legend=c("Female", "Male"), fill = c(my_red, my_blue), title="Пол")

На графиках видно, что со временем растёт количество писательниц.

2. Анализ данных с помощью cdplot, dotchart, boxplot, stripchart¶

cdplot¶

Для визуализации связи между двумя переменными, одна из которых является количественной, а другая качественной, можно использовать график условной плотности. Рассмотрим, например, какова связь между языком и количеством проданных копий.


Так как <<оценки условных плотностей более надёжны для областей с высокой плотностью $x$>>, а в данных наблюдается <<перекос>> (см. stripchart), то отображать на графике будем только часть данных - исключим выбросы справа (см. boxplot). В противном случае на части, где данных мало, график начинает выглядеть неадекватно.

Python:

In [15]:
lang_colors = reversed(["#B81840", "#C43142", "#D04544", "#DB5745", "#E56946",
                        "#EE7D47", "#F2994E", "#F5B355", "#F3A469", "#EC867E",
                        "#DD708E", "#BD6C99", "#9A68A4", "#76669B", "#4D6291", "#26456E"])
In [16]:
plt.figure(figsize=(10, 6))

ax = sns.kdeplot(data=df1[df1.sales < 66], x='sales', hue='language',
                 multiple='fill',
                 palette=lang_colors)

plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('Распределение количества проданных экземпляров по языкам', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.9), borderaxespad=0,
                title='Язык', title_fontsize=10, fontsize=10)
ax.set_xlim(min(df1.sales), 66)
plt.show()

R:

In [17]:
%%R

lang_colors <- c("#B81840", "#C43142", "#D04544", "#DB5745", "#E56946",
                 "#EE7D47", "#F2994E", "#F5B355", "#F3A469", "#EC867E",
                 "#DD708E", "#BD6C99", "#9A68A4", "#76669B", "#4D6291", "#26456E")
In [18]:
%%R -w 1000 -h 600

# добавляем дополнительное место слева от графика для легенды
par(mar=c(5, 4, 4, 16), xpd=TRUE)

languages <- factor(df1[df1$sales < 66, ]$language)
lang_levels <- unique(languages)
sales <- df1[df1$sales < 66, ]$sales

cdplot(languages ~ sales, col = lang_colors, ylevels = lang_levels,
       xlab = "Количество проденных экземпляров", ylab="", yaxlab = "",
       main = "Распределение количества проданных экземпляров по языкам")

# добавляем метки осей слева
#axis(2)

# добавляем легенду вне графика
legend("topright", inset=c(-0.2, 0.25), legend=lang_levels, fill = rev(lang_colors), title="Язык")

На графиках хорошо видно, что книги на английском языке чаще становятся бестселлерами.

dotchart¶

Одна из альтернатив столбчатым диаграммам - точечные диаграммы Кливленда. С их помощью рассмотрим, как связаны продажи книги с языком, на котором она написана. Для каждого языка найдём книгу с максимальным количеством проданных экземпляров и нанесём соответствующее значение на график (причём отсортируем языки в порядке возрастания количества произведений, написанных на том же языке).

Python:

In [19]:
language_freq = df1['language'].value_counts().reset_index()
language_freq.columns = ['language', 'freq']

language_max_sales = df1.groupby('language')['sales'].max().reset_index()
language_max_sales.columns = ['language', 'max_sales']

tmp1_df = pd.merge(language_freq, language_max_sales)
tmp1_df = tmp1_df.sort_values(['freq', 'max_sales'], ascending = (True, True))
In [20]:
plt.figure(figsize=(10, 6))

ax = sns.stripplot(x=tmp1_df.max_sales, y=tmp1_df.language, hue=tmp1_df.freq, palette='deep', size=7)

plt.xticks(fontsize=9)
plt.xlabel('Максимальное количество проданных экземпляров, млн', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Язык', fontsize=10)
plt.title('Связь максимального количества проданных экземпляров с количеством произведений на языке', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0, title='Количество произведений', fontsize=10, title_fontsize=10)
plt.show()

R:

In [21]:
%%R

language_freq <- as.data.frame(table(df1$language))
names(language_freq) <- c("language", "freq")

language_max_sales <- df1[order(df1$language, - df1$sales), ]
language_max_sales <- language_max_sales[!duplicated(language_max_sales$language), ]
language_max_sales <- language_max_sales[ , c("language", "sales")]

tmp1_df <- merge(language_max_sales, language_freq, by="language")
tmp1_df <- tmp1_df[order(-tmp1_df$freq, -tmp1_df$sales), ]
names(tmp1_df) <- c("language", "max_sales", "freq")

tmp1_df[tmp1_df$freq == 1, "color"] <- "royalblue"
tmp1_df[tmp1_df$freq == 2, "color"] <- "yellow"
tmp1_df[tmp1_df$freq == 3, "color"] <- "green"
tmp1_df[tmp1_df$freq == 4, "color"] <- "red"
tmp1_df[tmp1_df$freq == 5, "color"] <- "purple"
tmp1_df[tmp1_df$freq == 6, "color"] <- "brown"
tmp1_df[tmp1_df$freq == 132, "color"] <- "pink"
In [22]:
%%R -w 1000 -h 600

dotchart(tmp1_df$max_sales,
    labels = tmp1_df$language, groups = as.factor(tmp1_df$freq),
    pch = 16, cex = 1.2,
    color = tmp1_df$color, gcolor = "black",
    xlab = "Максимальное количество проданных экземпляров, млн", ylab = "Язык",
    main = "Связь максимального количества проданных экземпляров с количеством произведений на языке")

Видно, что с увеличением количества произведений, написанных на каком-либо языке, растёт и максимальные количество проданных экземпляров.

boxplot¶

Для того, чтобы визуализировать статистическую характеристику анализируемой совокупности (медиану, нижний и верхний квартили, максимум, минимум, выбросы), используются диаграммы размахов. Например, рассмотрим, как распределено количество проданных экземпляров.

Rython:

In [23]:
plt.figure(figsize=(10, 3))

sns.boxplot(data=df1, x='sales', width=0.4)

plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.title('Распределение количества проданных экземпляров', fontsize=10, fontweight='bold')
plt.show()

R:

In [24]:
%%R -w 1000 -h 300

boxplot(df1$sales, horizontal = TRUE,
        col = "royalblue",
        xlab = "Количество проданных экземпляров, млн",
        main = "Распределение количества проданных экземпляров")

С помощью диаграмм размахов можно сравнивать распределения разных величин. Для примера рассмотрим, чем отличаются распределения для разных языков.

Python:

In [25]:
tmp2_df = pd.merge(df1, tmp1_df)
tmp2_df = tmp2_df.sort_values(['freq', 'max_sales'], ascending = (True, True))
In [26]:
plt.figure(figsize=(10, 6))

ax = sns.boxplot(data=tmp2_df, x='sales', y='language', dodge=False, hue='freq', palette='deep')

plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Язык', fontsize=10)
plt.title('Распределение количества проданных экземпляров по языкам', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
                title='Количество произведений', fontsize=10, title_fontsize=10)
plt.show()

R:

In [27]:
%%R

deep_palette = rep(c("pink", "brown", "purple", "red", "green", "yellow", "royalblue"), times=c(1, 1, 3, 2, 1, 3, 5))
short_deep_palette = c("royalblue", "yellow", "green", "red", "purple", "brown", "pink")
language_levels = c('English', 'Russian', 'French', 'German', 'Japanese', 'Chinese', 'Italian', 'Spanish', 'Hindi', 'Norwegian', 'Swedish', 'Portuguese', 'Dutch', 'Czech', 'Gujarati', 'Yiddish')
In [28]:
%%R

tmp2_df <- data.frame(df1)
tmp2_df$language <- factor(tmp2_df$language, levels = language_levels)
In [29]:
%%R -w 1000 -h 600

# разворачиваем  метки осей горизонтально и добавляем дополнительное место справа от графика для легенды и слева для названия оси y
par(las=1, mar=c(5, 7, 4, 16), xpd=TRUE)

boxplot(tmp2_df$sales ~ tmp2_df$language,
        horizontal = TRUE, col=deep_palette,
        xlab = "Количество проданных экземпляров, млн",
        ylab="", # убираем название оси
        main = "Распределение количества проданных экземпляров")

# добавляем легенду вне графика
legend("topright", inset=c(-0.25, 0.35), legend=c(1,2,3,4,5,6,132), fill = short_deep_palette, title="Количество произведений")

# устанавливаем назание оси вручную, чтобы оно не накладывалось на метки оси
title(ylab = "Язык", line = 6)

На этих графиках зависимость количества проданных экземпляров от количества книг на языке не настолько явная, как казалось, когда рассматривали только максимумы (см. dotchart).

stripchart¶

Для визуализации распределения количественных переменных можно пользоваться одномерными диаграммами рассеивания. Рассмотрим, как распределено количество проданных копий.

Python:

In [30]:
plt.figure(figsize=(10, 3))

sns.stripplot(data=df1, x='sales', orient='h', jitter=True, alpha=0.6, size=7)

plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.xticks(fontsize=9)
plt.title('Распределение количества проданных экземпляров', fontsize=10, fontweight='bold')
plt.show()

R:

In [31]:
%%R -w 1000 -h 300

stripchart(df1$sales,
       xlab = "Количество проданных экземпляров, млн",
       vertical = FALSE,
       method = "jitter",
       pch = 19, col = my_blue, cex = 1.5,
       main = "Распределение количества проданных экземпляров")

С помощью диаграмм рассеивания можно сравнивать распределения разных величин. Для примера рассмортим, чем отличаются распределения для разных языков.

Python:

In [32]:
plt.figure(figsize=(10, 6))

ax = sns.stripplot(data=tmp2_df.query("language in ['English', 'Russian', 'French', 'German', 'Japanese', 'Chinese', 'Italian', 'Spanish', 'Hindi', 'Norwegian', 'Swedish', 'Portuguese', 'Dutch', 'Czech', 'Yiddish', 'Gujarati']"),
                      x='sales', y='language',
                      orient='h', dodge=False, jitter=True, alpha=0.6, size=7, hue='freq', palette='deep')

plt.xticks(fontsize=9)
plt.xlabel('Количество проданных экземпляров, млн', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Язык', fontsize=10)
plt.title('Распределение количества проданных экземпляров по языкам', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
                title='Количество произведений', fontsize=10, title_fontsize=10)
plt.show()

R:

In [33]:
%%R -w 1000 -h 600

# разворачиваем  метки осей горизонтально и добавляем дополнительное место справа от графика для легенды и слева для названия оси y
par(las=1, mar=c(5, 7, 4, 16), xpd=TRUE)

for(i in 1:length(deep_palette)){
    deep_palette[i] <- t_col(deep_palette[i], perc = 60)
}

for(i in 1:length(short_deep_palette)){
    short_deep_palette[i] <- t_col(short_deep_palette[i], perc = 60)
}

stripchart(tmp2_df$sales ~ tmp2_df$language,
           vertical = FALSE,
           method = "jitter",
           pch = 19, col = deep_palette, cex = 1.5,
           xlab = "Количество проданных экземпляров, млн",
           ylab="", # убираем название оси
           main = "Распределение количества проданных экземпляров")

# добавляем легенду вне графика
legend("topright", inset=c(-0.25, 0.35), legend=c(1,2,3,4,5,6,132), pch=19, col = short_deep_palette,, title="Количество произведений")

# устанавливаем назание оси вручную, чтобы оно не накладывалось на метки оси
title(ylab = "Язык", line = 6)

На графиках опять наблюдается рост количества проданных экземпляров с ростом количества книг на языке.

3. Проверка выбросов с точки зрения формальных статистических критериев Граббса и Q-теста Диксона¶

Для проверки выбросов с помощью этих критериев формулируются 2 гипотезы (нулевая и альтернативная):

$H_0:$ Подозрительное значение не является выбросом.

$H_1:$ Подозрительное значение - выброс.

Если полученый $pvalue$ меньше заданного заранее уровня значимости $\alpha$ (возьмём $\alpha = 0.05$), то $H_0$ должна быть отвергнута в пользу $H_1$, иначе - нет оснований отвергнуть нулевую гипотезу.

Также и критерий Граббса, и Q-тест Диксона основаны на предположении о нормальности распределения. Поэтому применим их к данным, близким к нормальным (см. график ядерных оценок плотности) - распределению годов первой публикации.

Критерий Граббса¶

Python:


В реализации теста Граббса на Python $pvalue$ "спрятан", поэтому выводим только результат.

In [34]:
min_out_index = grubbs.min_test_indices(df1.first_published, alpha=0.05)
In [35]:
max_out_index = grubbs.max_test_indices(df1.first_published, alpha=0.05)
In [36]:
df1.iloc[min_out_index + max_out_index]
Out[36]:
title author gender birth_year age language first_published sales genres
31 The Divine Comedy (La Divina Commedia) Dante Alighieri M 1265 39 Italian 1304 11.5 [Classics, Poetry, Fiction, Literature, Philos...

R:

In [37]:
%%R

res <- grubbs.test(df1[]$first_published, type = 10)

print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет")
	Grubbs test for one outlier

data:  df1[]$first_published
G = 10.27459, U = 0.38981, p-value < 2.2e-16
alternative hypothesis: lowest value 1304 is an outlier

Есть ли основания отбросить гипотезу? Да
In [38]:
%%R

res <- grubbs.test(df1[]$first_published, type = 10, opposite = TRUE)

print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет")
	Grubbs test for one outlier

data:  df1[]$first_published
G = 0.86350, U = 0.99569, p-value = 1
alternative hypothesis: highest value 2018 is an outlier

Есть ли основания отбросить гипотезу? Нет
In [39]:
%%R

df1[df1$first_published == 1304, ]
                                    title          author gender birth_year age
32 The Divine Comedy (La Divina Commedia) Dante Alighieri      M       1265  39
   language first_published sales
32  Italian            1304  11.5
                                                                 genres
32 Classics, Poetry, Fiction, Literature, Philosophy, Religion, Fantasy

Тест Граббса нашёл один выброс - "Божественную комедию" Данте Алигьери (1304).

Q-тест Диксона¶

Тест Диксона может обрабатывать выборки размером не более 30. Поэтому выделим диапазон из не более чем 30 значений, в который потенциально может попасть выброс. Также если в диапазон не попала часть одинаковых значений, удалим все такие значения.

Python:

In [40]:
min_df = df1.sort_values(['first_published']).reset_index()[0:31]
min_df = min_df[min_df.first_published != min_df.first_published[30]]
print('Количество элементов для проверки минимума:', min_df.shape[0])
Количество элементов для проверки минимума: 30
In [41]:
max_df = df1.sort_values(['first_published'], ascending=False).reset_index()[0:31]
max_df = max_df[max_df.first_published != max_df.first_published[30]]
print('Количество элементов для проверки максимума:', max_df.shape[0])
Количество элементов для проверки максимума: 29

R:

In [42]:
%%R

min_df <- df1[order(df1$first_published, decreasing=FALSE), ]$first_published[1:31]
min_df <- min_df[min_df != min_df[31]]
cat("Количество элементов для проверки минимума:", length(min_df))
Количество элементов для проверки минимума: 30
In [43]:
%%R

max_df <- df1[order(df1$first_published, decreasing=TRUE), ]$first_published[1:31]
max_df <- max_df[max_df != max_df[31]]
cat("Количество элементов для проверки максимума:", length(max_df))
Количество элементов для проверки максимума: 29

Для минимальных и максимальных выбросов получились выборки размером 30 и 29 соответственно.

Python:


На Python нет готовой реализации теста Диксона, поэтому сделаем её сами. Для выборок размера 30 и 29 при уровне значимости 0.95 критические значения равны 0.290 и 0.301 соответственно (это табличные данные, полученные моделированием; численно получить их можно только для выборки размера $n=3$).

In [44]:
out = []

maX = min_df.first_published[29]
miN = min_df.first_published[0]
second_min = min(min_df[min_df.first_published != miN].first_published)
min_q = (second_min - miN) / (maX - miN)
if (min_q >= 0.290):
    out.append(miN)

maX = max_df.first_published[0]
second_max = max(max_df[max_df.first_published != maX].first_published)
miN = max_df.first_published[28]
max_q = (maX - second_max) / (maX - miN)
if (max_q >= 0.310):
    out.append(maX)

df1[df1.first_published.isin(out)]
Out[44]:
title author gender birth_year age language first_published sales genres
31 The Divine Comedy (La Divina Commedia) Dante Alighieri M 1265 39 Italian 1304 11.5 [Classics, Poetry, Fiction, Literature, Philos...

R:

In [45]:
%%R

res <- dixon.test(min_df)

print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет")
	Dixon test for outliers

data:  min_df
Q = 0.77057, p-value < 2.2e-16
alternative hypothesis: lowest value 1304 is an outlier

Есть ли основания отбросить гипотезу? Да
In [46]:
%%R

res <- dixon.test(max_df)

print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет")
	Dixon test for outliers

data:  max_df
Q = 0.13333, p-value = 0.6394
alternative hypothesis: highest value 2018 is an outlier

Есть ли основания отбросить гипотезу? Нет

Тест Диксона дал тот же результат, что и тест Граббса. Визуализируем его.

Python:

In [47]:
plt.figure(figsize=(10, 3))

ax = sns.stripplot(data=df1[df1.first_published.isin(out)], x='first_published', orient='h', jitter=True,
              alpha=0.6, size=7, color='red', label='Выброс')
sns.stripplot(data=df1[~df1.first_published.isin(out)], x='first_published', orient='h', jitter=True,
              alpha=0.6, size=7, color='green', label='Остальная выборка')

plt.xticks(fontsize=9)
plt.xlabel('Год первой публикации', fontsize=10)
plt.title('Распределение первых публикаций', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
                title='')
plt.show()

R:

In [48]:
%%R -w 1000 -h 300

my_green <- t_col("forestgreen", perc = 60)
stripchart(df1[df1$first_published != 1304, ]$first_published,
       method = "jitter",
       vertical = FALSE,
       pch = 19, col = my_green, cex = 1.5,
       xlim = c(1304, 2018),
       xlab = "Распределение первых публикаций",
       main = "Распределение количества проданных экземпляров")

stripchart(df1[df1$first_published == 1304, ]$first_published,
       method = "jitter",
       pch = 19, col = my_red, cex = 1.5, add = TRUE)

4. Заполнение пропусков в данных¶

Рассмотрим распределение возрастов авторов.

Python:

In [49]:
plt.figure(figsize=(10, 3))

sns.stripplot(data=df1, x='age', orient='h', jitter=True,
              alpha=0.6, size=7)

plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.title('Распределение возрастов авторов', fontsize=10, fontweight='bold')
plt.show()

R:

In [50]:
%%R -w 1000 -h 300

stripchart(df1$age,
       xlab = "Возраст",
       vertical = FALSE,
       method = "jitter",
       pch = 19, col = my_blue, cex = 1.5,
       main = "Распределение возрастов авторов")

Удалим возраст в некоторых строках.

Python:

In [51]:
df_copy = df1.sort_values(['age', 'author']).reset_index()
df_copy.iloc[[36, 89, 141, 147]]
Out[51]:
index title author gender birth_year age language first_published sales genres
36 104 Anne of Green Gables Lucy Maud Montgomery F 1874 34 English 1908 50.0 [Classics, Fiction, Young Adult, Historical Fi...
89 165 Man's Search for Meaning (Ein Psychologe erleb... Viktor Frankl M 1905 41 German 1946 12.0 [Nonfiction, Psychology, Philosophy, History, ...
141 18 The Lion, the Witch and the Wardrobe C. S. Lewis M 1898 52 English 1950 85.0 [Fantasy, Classics, Fiction, Young Adult, Chil...
147 85 Heidi Johanna Spyri F 1827 53 German 1880 50.0 [Classics, Fiction, Childrens, Young Adult, Hi...
In [52]:
df_copy['age'][36] = np.nan
df_copy['age'][89] = np.nan
df_copy['age'][141] = np.nan
df_copy['age'][147] = np.nan

df_copy['age'].isna().sum()
Out[52]:
4

R:

In [53]:
%%R

df_copy <- df1[order(c(df1$age)), ]
df_copy[c(37, 90, 142, 148), ]
                                                                       title
105                                                     Anne of Green Gables
166 Man's Search for Meaning (Ein Psychologe erlebt das Konzentrationslager)
19                                      The Lion, the Witch and the Wardrobe
86                                                                     Heidi
                  author gender birth_year age language first_published sales
105 Lucy Maud Montgomery      F       1874  34  English            1908    50
166        Viktor Frankl      M       1905  41   German            1946    12
19           C. S. Lewis      M       1898  52  English            1950    85
86         Johanna Spyri      F       1827  53   German            1880    50
                                                                                  genres
105  Classics, Fiction, YoungAdult, HistoricalFiction, Childrens, MiddleGrade, Audiobook
166             Nonfiction, Psychology, Philosophy, History, SelfHelp, Memoir, Biography
19             Fantasy, Classics, Fiction, YoungAdult, Childrens, MiddleGrade, Adventure
86  Classics, Fiction, Childrens, YoungAdult, HistoricalFiction, MiddleGrade, Historical
In [54]:
%%R

df_copy$age[37] <- NA
df_copy$age[90] <- NA
df_copy$age[142] <- NA
df_copy$age[148] <- NA

sum(is.na(df_copy$age))
[1] 4

Заполним пропуски средним по всем значениям и cравним получившееся распределение с истинными.

Python:

In [55]:
# одинаковые границы отрезков разбиения для гистограмм
edges = np.histogram_bin_edges(df1['age'], bins=14, range=(15,85))
edges
Out[55]:
array([15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75.,
       80., 85.])
In [56]:
all_mean_df = df_copy.copy()

all_mean = all_mean_df.mean()['age']

all_mean_df['age'][36] = all_mean
all_mean_df['age'][89] = all_mean
all_mean_df['age'][141] = all_mean
all_mean_df['age'][147] = all_mean

print(all_mean)
42.7953216374269
In [57]:
fig = plt.figure(figsize=(10,3))

ax1 = plt.subplot(121)
sns.histplot(data=df1.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')

ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
sns.histplot(data=all_mean_df.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены средним по всем значениям', fontsize=10, fontweight='bold')

plt.show()

R:

In [58]:
%%R

# одинаковые границы отрезков разбиения для гистограмм
edges = hist(df1$age, freq = FALSE, breaks = "FD", plot = FALSE)
edges[1]
$breaks
 [1] 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85

In [59]:
%%R

all_mean_df <- data.frame(df_copy)

all_mean <- mean(all_mean_df$age, na.rm = TRUE)

all_mean_df$age[37] <- all_mean
all_mean_df$age[90] <- all_mean
all_mean_df$age[142] <- all_mean
all_mean_df$age[148] <- all_mean

cat(all_mean)
42.79532
In [60]:
%%R -w 1000 -h 300

par(mfrow = c(1,2))

hist(df1$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Истинные значения")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)

hist(all_mean_df$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены средним по всем значениям")
lines(density(all_mean_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)

Заполним пропуски медианой по всем значениям и вычислим отклонение заполнений от истинных значений.

Python:

In [61]:
all_median_df = df_copy.copy()

all_median = all_median_df.median()['age']

all_median_df['age'][36] = all_median
all_median_df['age'][89] = all_median
all_median_df['age'][141] = all_median
all_median_df['age'][147] = all_median

print(all_median)
41.0
In [62]:
fig = plt.figure(figsize=(10,3))

ax1 = plt.subplot(121)
sns.histplot(data=df1.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')

ax2 = plt.subplot(122, sharex=ax1, sharey=ax1)
sns.histplot(data=all_median_df.age, kde=True, bins=edges, alpha=0.3, stat='density')
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены медианой по всем значениям', fontsize=10, fontweight='bold')

plt.show()

R:

In [63]:
%%R

all_median_df <- data.frame(df_copy)

all_median <- median(all_median_df$age, na.rm = TRUE)

all_median_df$age[37] <- all_median
all_median_df$age[90] <- all_median
all_median_df$age[142] <- all_median
all_median_df$age[148] <- all_median

cat(all_median)
41
In [64]:
%%R -w 1000 -h 300

par(mfrow = c(1,2))

hist(df1$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Истинные значения")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)

hist(all_median_df$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены медианой по всем значениям")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)

Пропуски можно заполнять средними и медианами не по всем значениям, а по категории. Для примера рассмотрим, чем отличаются распределения для женщин и мужчин.

Python:

In [65]:
plt.figure(figsize=(10, 3))

ax = sns.stripplot(data=df1, x='age', orient='h', jitter=True, dodge=True, hue='gender',
              alpha=0.4, size=7, palette=['red', 'blue'])

plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.title('Распределение возрастов авторов', fontsize=10, fontweight='bold')
sns.move_legend(ax, loc='upper left', bbox_to_anchor=(1.05, 0.65), borderaxespad=0,
                title='Пол автора', fontsize=10, title_fontsize=10)
plt.show()

R:

In [66]:
%%R -w 1000 -h 300

# добавляем дополнительное место справа от графика для легенды
par(mar=c(5, 4, 4, 16), xpd=TRUE)

stripchart(df1$age ~ df1$gender,
           vertical = FALSE,
           method = "jitter",
           pch = 19, col = c(my_red, my_blue), cex = 1.5,
           xlab = "Возраст",
           ylab="", yaxt="n",
           main = "Распределение возрастов авторов")

# добавляем легенду вне графика
legend("topright", inset=c(-0.15, 0.35), legend=c("M", "F"), pch=19,
       col = c(my_blue, my_red), title = "Пол автора")

Заполним пропуски средними по категориям и вычислим отклонение заполнений от истинных значений.

Python:

In [67]:
gender_mean_df = df_copy.copy()

m_mean = gender_mean_df[gender_mean_df.gender == 'M'].mean()['age']
f_mean = gender_mean_df[gender_mean_df.gender == 'F'].mean()['age']

gender_mean_df['age'][36] = m_mean if (gender_mean_df['gender'][36] == 'M') else f_mean
gender_mean_df['age'][89] = m_mean if (gender_mean_df['gender'][89] == 'M') else f_mean
gender_mean_df['age'][141] = m_mean if (gender_mean_df['gender'][141] == 'M') else f_mean
gender_mean_df['age'][147] = m_mean if (gender_mean_df['gender'][147] == 'M') else f_mean

print(m_mean, f_mean)
43.4 41.55357142857143
In [68]:
fig = plt.figure(figsize=(10,9))

ax1 = plt.subplot(321)
sns.histplot(data=df1.age,
             kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')

ax2 = plt.subplot(322, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_mean_df.age,
             kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены средними по категориям', fontsize=10, fontweight='bold')

ax3 = plt.subplot(323, sharex=ax1, sharey=ax1)
sns.histplot(data=df1[df1.gender == 'M'].age,
             kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax3.set_title('Истинные значения (муж)', fontsize=10, fontweight='bold')

ax4 = plt.subplot(324, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_mean_df[gender_mean_df.gender == 'M'].age,
             kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax4.set_title('Пропуски заполнены средними по категориям (муж)', fontsize=10, fontweight='bold')

ax5 = plt.subplot(325, sharex=ax1, sharey=ax1)
sns.histplot(data=df1.age[df1.gender == 'F'],
             kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax5.set_title('Истинные значения (жен)', fontsize=10, fontweight='bold')

ax6 = plt.subplot(326, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_mean_df[gender_mean_df.gender == 'F'].age,
             kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax6.set_title('Пропуски заполнены средними по категориям (жен)', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()

R:

In [69]:
%%R

gender_mean_df <- data.frame(df_copy)

m_mean <- mean(gender_mean_df[gender_mean_df$gender %like% "M", "age"], na.rm = TRUE)
f_mean <- mean(gender_mean_df[gender_mean_df$gender %like% "F", "age"], na.rm = TRUE)

gender_mean_df$age[37] <- if (gender_mean_df$gender[37] == "M") m_mean else f_mean
gender_mean_df$age[90] <- if (gender_mean_df$gender[90] == "M") m_mean else f_mean
gender_mean_df$age[142] <- if (gender_mean_df$gender[142] == "M") m_mean else f_mean
gender_mean_df$age[148] <- if (gender_mean_df$gender[148] == "M") m_mean else f_mean

cat(m_mean, f_mean)
43.4 41.55357
In [70]:
%%R -w 1000 -h 900

par(mfrow = c(3,2))

hist(df1$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Истинные значения")
lines(density(df1$age, adjust=), col = "blue", lwd = 3)

hist(gender_mean_df$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены средними по категориям")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "blue", lwd = 3)

hist(df1[df1$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_blue,
     xlab = "Возраст", ylab = "",
     main = "Истинные значения (муж)")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)

hist(gender_mean_df[gender_mean_df$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_blue,
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены средними по категориям (муж)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)

hist(df1[df1$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_red,
     xlab = "Возраст", ylab = "",
     main = "Истинные значения (жен)")
lines(density(df1$age, adjust=), col = "red", lwd = 3)

hist(gender_mean_df[gender_mean_df$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_red,
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены средними по категориям (жен)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "red", lwd = 3)

Заполним пропуски медианами по категориям и вычислим отклонение заполнений от истинных значений.

Python:

In [71]:
gender_median_df = df_copy.copy()

m_median = gender_median_df[gender_median_df.gender == 'M'].median()['age']
f_median = gender_median_df[gender_median_df.gender == 'F'].median()['age']

gender_median_df['age'][36] = m_median if (gender_median_df['gender'][36] == 'M') else f_median
gender_median_df['age'][89] = m_median if (gender_median_df['gender'][89] == 'M') else f_median
gender_median_df['age'][141] = m_median if (gender_median_df['gender'][141] == 'M') else f_median
gender_median_df['age'][147] = m_median if (gender_median_df['gender'][147] == 'M') else f_median

print(m_median, f_median)
41.0 40.0
In [72]:
fig = plt.figure(figsize=(10,9))

ax1 = plt.subplot(321)
sns.histplot(data=df1.age,
             kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax1.set_title('Истинные значения', fontsize=10, fontweight='bold')

ax2 = plt.subplot(322, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_median_df.age,
             kde=True, bins=edges, alpha=0.3, stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax2.set_title('Пропуски заполнены медианами по категориям', fontsize=10, fontweight='bold')

ax3 = plt.subplot(323, sharex=ax1, sharey=ax1)
sns.histplot(data=df1[df1.gender == 'M'].age,
             kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax3.set_title('Истинные значения (муж)', fontsize=10, fontweight='bold')

ax4 = plt.subplot(324, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_median_df[gender_median_df.gender == 'M'].age,
             kde=True, bins=edges, alpha=0.3, color='blue', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax4.set_title('Пропуски заполнены медианами по категориям (муж)', fontsize=10, fontweight='bold')

ax5 = plt.subplot(325, sharex=ax1, sharey=ax1)
sns.histplot(data=df1.age[df1.gender == 'F'],
             kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax5.set_title('Истинные значения (жен)', fontsize=10, fontweight='bold')

ax6 = plt.subplot(326, sharex=ax1, sharey=ax1)
sns.histplot(data=gender_median_df[gender_median_df.gender == 'F'].age,
             kde=True, bins=edges, alpha=0.3, color='red', stat='density')
plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
ax6.set_title('Пропуски заполнены медианами по категориям (жен)', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()

R:

In [73]:
%%R

gender_median_df <- data.frame(df_copy)

m_median <- median(gender_median_df[gender_median_df$gender %like% "M", "age"], na.rm = TRUE)
f_median <- median(gender_median_df[gender_median_df$gender %like% "F", "age"], na.rm = TRUE)

gender_median_df$age[37] <- if (gender_median_df$gender[37] == "M") m_median else f_median
gender_median_df$age[90] <- if (gender_median_df$gender[90] == "M") m_median else f_median
gender_median_df$age[142] <- if (gender_median_df$gender[142] == "M") m_median else f_median
gender_median_df$age[148] <- if (gender_median_df$gender[148] == "M") m_median else f_median

cat(m_median, f_median)
41 40
In [74]:
%%R -w 1000 -h 900

par(mfrow = c(3,2))

hist(df1$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Истинные значения")
lines(density(df1$age, adjust=), col = "blue", lwd = 3)

hist(gender_median_df$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены медианами по категориям")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "blue", lwd = 3)

hist(df1[df1$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_blue,
     xlab = "Возраст", ylab = "",
     main = "Истинные значения (муж)")
lines(density(df1$age, adjust=), col = "royalblue", lwd = 3)

hist(gender_median_df[gender_median_df$gender %like% "M", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_blue,
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены медианами по категориям (муж)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "royalblue", lwd = 3)

hist(df1[df1$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_red,
     xlab = "Возраст", ylab = "",
     main = "Истинные значения (жен)")
lines(density(df1$age, adjust=), col = "red", lwd = 3)

hist(gender_median_df[gender_median_df$gender %like% "F", ]$age, freq = FALSE, breaks = edges$breaks,
     col = my_red,
     xlab = "Возраст", ylab = "",
     main = "Пропуски заполнены медианами по категориям (жен)")
lines(density(all_median_df$age, na.rm = TRUE, adjust=), col = "red", lwd = 3)

Сравним разные способы заполнения.

Python:

In [75]:
real_age = df1.sort_values(['age']).reset_index()[['age']].iloc[[36, 89, 141, 147]]

all_mean_age = all_mean_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]
all_median_age = all_median_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]
gender_mean_age = gender_mean_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]
gender_median_age = gender_median_df.reset_index()[['age']].iloc[[36, 89, 141, 147]]

all_mean_age['age'] = real_age['age'] - all_mean_age['age']
all_median_age['age'] = real_age['age'] - all_median_age['age']
gender_mean_age['age'] = real_age['age'] - gender_mean_age['age']
gender_median_age['age'] = real_age['age'] - gender_median_age['age']
In [76]:
fig = plt.figure(figsize=(10,6))

ax1 = plt.subplot(221)
plt.bar(data=all_mean_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
        height=all_mean_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены средним по всем', fontsize=10, fontweight='bold')

ax2 = plt.subplot(222, sharex=ax1, sharey=ax1)
plt.bar(data=all_median_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
        height=all_median_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены медианой по всем', fontsize=10, fontweight='bold')

ax3 = plt.subplot(223, sharex=ax1, sharey=ax1)
plt.bar(data=gender_mean_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
        height=gender_mean_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены средними по категориям', fontsize=10, fontweight='bold')

ax4 = plt.subplot(224, sharex=ax1, sharey=ax1)
plt.bar(data=gender_median_age.query("index in ['36', '89', '141', '147']"), x=['36', '89', '141', '147'],
        height=gender_median_age.age)
plt.xticks(fontsize=9)
plt.xlabel('Индекс пропуска', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Разница с истинным значением', fontsize=10)
plt.title('Пропуски заполнены медианами по категориям', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=3)

plt.show()
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.
Using categorical units to plot a list of strings that are all parsable as floats or dates. If these strings should be plotted as numbers, cast to the appropriate data type before plotting.

R:

In [77]:
%%R

real_age <- df1[order(c(df1$age)), ]$age[c(37, 90, 142, 148)]

all_mean_age <- all_mean_df$age[c(37, 90, 142, 148)]
all_median_age <- all_median_df$age[c(37, 90, 142, 148)]
gender_mean_age <- gender_mean_df$age[c(37, 90, 142, 148)]
gender_median_age <- gender_median_df$age[c(37, 90, 142, 148)]

all_mean_age <- real_age - all_mean_age
all_median_age<- real_age - all_median_age
gender_mean_age <- real_age - gender_mean_age
gender_median_age <- real_age - gender_median_age
In [78]:
%%R -w 1000 -h 600

par(mfrow = c(2,2))

barplot(all_mean_age, names.arg = c("37", "90", "142", "148"),
     col = "lightblue",
     xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
     main = "Пропуски заполнены средним по всем")


barplot(all_median_age, names.arg = c("37", "90", "142", "148"),
     col = "lightblue",
     xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
     main = "Пропуски заполнены медианой по всем")

barplot(gender_mean_age, names.arg = c("37", "90", "142", "148"),
     col = "lightblue",
     xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
     main = "Пропуски заполнены средними по категориям")

barplot(gender_median_age, names.arg = c("37", "90", "142", "148"),
     col = "lightblue",
     xlab = "Индекс пропуска", ylab = "Разница с истинным значением",
     main = "Пропуски заполнены медианами по категориям")

В данном случае разные способы заполнения пропусков дают примерно одинаковый результат, но при заполнениях медианами значения получаются меньше из-за ассиметрии данных.

5. Анализ данных с помощью графиков квантилей, метода огибающих и стандартных процедур проверки гипотез о нормальности (сгенерированные данные из нормального распределения)¶

Сгенерируем данные из нормальных распределений с различными параметрами:

  • Выборки малого размера (100 наблюдений):
    • N(0, 1)
    • N(0, 2)
  • Выборки среднего размера (1000 наблюдений):
    • N(0, 1)
    • N(2, 1)

Python:

In [79]:
data1 = np.random.normal(0, 1, size=100)
data2 = np.random.normal(0, 2, size=100)
data3 = np.random.normal(0, 1, size=1000)
data4 = np.random.normal(2, 1, size=1000)

R:

In [80]:
%%R

data1 <- rnorm(100, mean = 0, sd = 1)
data2 <- rnorm(100, mean = 0, sd = 2)
data3 <- rnorm(1000, mean = 0, sd = 1)
data4 <- rnorm(1000, mean = 2, sd = 1)

Эмпирические функции распределения¶

Для проверки нормальности можно использовать эмпирические функции распределения, сравнивая их графики с графиками теоретических функций распределения. В случае нормального распределения они должны быть похожи.

Python:

In [81]:
fig = plt.figure(figsize=(10,6))
x = np.linspace(-10, 10, 1000)

ax1 = plt.subplot(221)
plt.plot(x, stats.norm.cdf(x), color='red')
sns.ecdfplot(data=data1)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax1.set_title('N(0,1), 100 наблюдений', fontsize=10, fontweight='bold')

ax2 = plt.subplot(222, sharex=ax1, sharey=ax1)
plt.plot(x, stats.norm.cdf(x, 0, 2), color='red')
sns.ecdfplot(data=data2)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax2.set_title('N(0,2), 100 наблюдений', fontsize=10, fontweight='bold')

ax3 = plt.subplot(223, sharex=ax1, sharey=ax1)
plt.plot(x, stats.norm.cdf(x), color='red')
sns.ecdfplot(data=data3)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax3.set_title('N(0,1), 1000 наблюдений', fontsize=10, fontweight='bold')

ax4 = plt.subplot(224, sharex=ax1, sharey=ax1)
plt.plot(x, stats.norm.cdf(x, 2, 1), color='red')
sns.ecdfplot(data=data4)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax4.set_title('N(2,1), 1000 наблюдений', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()

R:

In [82]:
%%R -w 1000 -h 600

par(mfrow = c(2,2))

plot(sort(data1), pnorm(sort(data1), mean = 0, sd = 1), type = "l",
     xlim = c(-10, 10), ylim = c(-0.01, 1.01),
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="N(0,1), 100 наблюдений")
plot(ecdf(data1), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

plot(sort(data2), pnorm(sort(data2), mean = 0, sd = 2), type = "l",
     xlim = c(-10, 10), ylim = c(-0.01, 1.01),
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="N(0,2), 100 наблюдений")
plot(ecdf(data2), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

plot(sort(data3), pnorm(sort(data3), mean = 0, sd = 1), type = "l",
     xlim = c(-10, 10), ylim = c(-0.01, 1.01),
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="N(0,1), 1000 наблюдений")
plot(ecdf(data3), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

plot(sort(data4), pnorm(sort(data4), mean = 2, sd = 1), type = "l",
     xlim = c(-10, 10), ylim = c(-0.01, 1.01),
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="N(2,1), 1000 наблюдений")
plot(ecdf(data4), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

Квантили¶

Проверять, получены ли данные из некоторого теоретического распределения, можно с помощью графиков QQ ("квантиль-квантиль"). В случае нормальных распределений с различными параметрами точки должны лежать вдоль соответствующих прямых.

Python:

In [83]:
fig = plt.figure(figsize=(10,6))

ax1 = plt.subplot(221)
sm.qqplot(data1, ax=ax1, line='45')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
ax1.set_xlim(-6, 6)
plt.yticks(fontsize=9)
ax1.set_ylim(-6, 6)
plt.ylabel('Выборочные квантили', fontsize=10)
ax1.set_title('N(0,1), 100 наблюдений', fontsize=10, fontweight='bold')

ax2 = plt.subplot(222, sharex=ax1, sharey=ax1)
sm.qqplot(data2, ax=ax2, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax2.set_title('N(0,2), 100 наблюдений', fontsize=10, fontweight='bold')

ax3 = plt.subplot(223, sharex=ax1, sharey=ax1)
sm.qqplot(data3, ax=ax3, line='45')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax3.set_title('N(0,1), 1000 наблюдений', fontsize=10, fontweight='bold')

ax4 = plt.subplot(224, sharex=ax1, sharey=ax1)
sm.qqplot(data4, ax=ax4, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax4.set_title('N(2,1), 1000 наблюдений', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()

R:

In [84]:
%%R -w 1000 -h 600

par(mfrow = c(2,2))

qqnorm(data1,
       xlim = c(-6, 6), ylim = c(-6, 6),
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "N(0,1), 100 наблюдений")
qqline(data1,
       col = "red", lwd = 2)

qqnorm(data2,
       xlim = c(-6, 6), ylim = c(-6, 6),
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "N(0,2), 100 наблюдений")
qqline(data2,
       col = "red", lwd = 2)

qqnorm(data3,
       xlim = c(-6, 6), ylim = c(-6, 6),
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "N(0,1), 1000 наблюдений")
qqline(data3,
       col = "red", lwd = 2)

qqnorm(data4,
       xlim = c(-6, 6), ylim = c(-6, 6),
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "N(2,1), 1000 наблюдений")
qqline(data4,
       col = "red", lwd = 2)

Метод огибающих¶

Ещё один способ визуально проверить нормальность распределения - метод огибающих.

Python:

In [85]:
def envelmet(x, title, ax):
    # точки
    z = (x - x.mean()) / x.std()
    x_qq, _ = stats.probplot(z)
    ax.plot(x_qq[0], x_qq[1], '.k')

    # линия
    ax.plot((-3, 3), (-3, 3), c='red')

    ax.tick_params(axis='both', labelsize=9)
    ax.set_xlabel('Квантили', fontsize=10)
    ax.set_ylabel('Статистики', fontsize=10)
    ax.set_title(title, fontsize=10, fontweight='bold')

fig, ax = plt.subplots(nrows=2, ncols=2, figsize = (10, 6))

envelmet(data1, 'N(0,1), 100 наблюдений', ax[0, 0])
envelmet(data2, 'N(0,2), 100 наблюдений', ax[0, 1])
envelmet(data3, 'N(0,1), 1000 наблюдений', ax[1, 0])
envelmet(data4, 'N(2,1), 1000 наблюдений', ax[1, 1])


fig.tight_layout()

plt.show()

R:

In [86]:
%%R -w 1000 -h 600

par(mfrow = c(2,2))

envelmet <- function(x, title){
    # точки
    z <- (x - mean(x)) / sqrt(var(x))
    x_qq <- qqnorm(z, plot.it = FALSE)
    x_qq <- lapply(x_qq, sort)
    plot(x_qq,
         xlim = c(-3, 3), ylim = c(-3, 3),
         xlab = "Квантили", ylab = "Статистики",
         main = title)

    # огибающие
    x.gen <- function(dat, mle) rnorm(length(dat))
    x.qqboot <- boot(z, sort, R = 999,
        sim = "parametric", ran.gen = x.gen)
    sapply(1:999, function(i) lines(x_qq$x, x.qqboot$t[i,], type = "l", col = "grey"))
    points (x_qq, pch = 20)

    #прямая
    lines(c(-3, 3), c(-3, 3), col = "red", lwd = 2)
}

envelmet(data1, "N(0,1), 100 наблюдений")
envelmet(data2, "N(0,2), 100 наблюдений")
envelmet(data3, "N(0,1), 1000 наблюдений")
envelmet(data4, "N(2,1), 1000 наблюдений")

Стандартные процедуры проверки нормальности¶

Для каждого из критериев формулируется нулевая гипотеза $H_0$ о характере распределения выборки $X_1, \dots, X_n$. Если полученный при использовании критерия $pvalue$ меньше заданного заранее уровня значимости $\alpha$ (положим $\alpha = 0.5$), то $H_0$ должна быть отвергнута, иначе - нет оснований отбросить нулевую гипотезу.

Критерий Колмогорова - Смирнова¶

$$H_0 : F_n(x) = F(x, \theta_0), \quad x \in \mathbb{R},$$

т.е.данные получены из распределения с заданными параметрами.

Положим $F(x, \theta_0) \equiv \Phi(x)$, чтобы проверить, является ли распределение стандартным нормальным.

Python:

In [87]:
_, pvalue = kstest(data1, 'norm')
print('N(0,1), 100 наблюдений:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')

_, pvalue = kstest(data2, 'norm')
print('N(0,2), 100 наблюдений:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')

_, pvalue = kstest(data3, 'norm')
print('N(0,1), 1000 наблюдений:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')


_, pvalue = kstest(data4, 'norm')
print('N(2,1), 1000 наблюдений:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
    pvalue = 0.48487951108184757
    Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
    pvalue = 0.0002731978404077222
    Есть ли основания отбросить гипотезу? Да
N(0,1), 1000 наблюдений:
    pvalue = 0.6034541166986684
    Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
    pvalue = 0.0
    Есть ли основания отбросить гипотезу? Да

R:

In [88]:
%%R

res <- ks.test(data1, "pnorm")
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- ks.test(data2, "pnorm")
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- ks.test(data3, "pnorm")
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- ks.test(data4, "pnorm")
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------
	Asymptotic one-sample Kolmogorov-Smirnov test

data:  data1
D = 0.071778, p-value = 0.6815
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,2), 100 наблюдений ------------------
	Asymptotic one-sample Kolmogorov-Smirnov test

data:  data2
D = 0.22247, p-value = 0.0001005
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Да 

------------------ N(0,1), 1000 наблюдений ------------------
	Asymptotic one-sample Kolmogorov-Smirnov test

data:  data3
D = 0.014548, p-value = 0.984
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Нет 

------------------ N(2,1), 1000 наблюдений ------------------
	Asymptotic one-sample Kolmogorov-Smirnov test

data:  data4
D = 0.69077, p-value < 2.2e-16
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Да 

Критерий Шапиро - Уилка¶

$$H_0: F_n(x) \in \{ \; F(x,\theta) \quad | \quad \theta \in \Theta_0 \; \} \equiv $$$$ \equiv \{ \; F(x; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \int^x_{-\infty} e^{-\frac{(t-\mu)^2}{2\sigma^2}} \; | \; \mu \in \mathbb{R}, \sigma > 0 \; \},$$

т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.

Также для критерия Шапиро-Уилка необходимо, чтобы размер выборки был не более 5000.

Python:

In [89]:
_, pvalue = shapiro(data1)
print('N(0,1), 100 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')

_, pvalue = shapiro(data2)
print('N(0,2), 100 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')

_, pvalue = shapiro(data3)
print('N(0,1), 1000 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')

_, pvalue = shapiro(data4)
print('N(2,1), 1000 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
    pvalue =  0.6277329921722412
    Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
    pvalue =  0.04325342923402786
    Есть ли основания отбросить гипотезу? Да 
N(0,1), 1000 наблюдений:
    pvalue =  0.9437151551246643
    Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
    pvalue =  0.2394869327545166
    Есть ли основания отбросить гипотезу? Нет

R:

In [90]:
%%R

res <- shapiro.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- shapiro.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- shapiro.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- shapiro.test(data4)
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------
	Shapiro-Wilk normality test

data:  data1
W = 0.99228, p-value = 0.841

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,2), 100 наблюдений ------------------
	Shapiro-Wilk normality test

data:  data2
W = 0.99281, p-value = 0.8764

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,1), 1000 наблюдений ------------------
	Shapiro-Wilk normality test

data:  data3
W = 0.99847, p-value = 0.5338

Есть ли основания отбросить гипотезу? Нет 

------------------ N(2,1), 1000 наблюдений ------------------
	Shapiro-Wilk normality test

data:  data4
W = 0.99719, p-value = 0.0786

Есть ли основания отбросить гипотезу? Нет 

Критерий Андерсона - Дарлинга¶

$$H_0: F_n(x) \in \{ \; F(x,\theta) \; | \; \theta \in \Theta_0 \; \} \equiv $$$$ \equiv \{ \; F(x; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \int^x_{-\infty} e^{-\frac{(t-\mu)^2}{2\sigma^2}} \; | \; \mu \in \mathbb{R}, \sigma > 0 \; \},$$

т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.

Python:


Так как в Python реализация критерия Андерсона-Дарлинга (scipy.stats.anderson()) не возвращает pvalue, нужно сравнить полученную тестовую статистику с табличным критическим значением (для уровня значимости $\alpha = 0.05$ и выборок размера 100 и 1000 критическое значение равно 0.759 и 0.784 соответственно).

In [91]:
res = anderson(data1, 'norm')
print('N(0,1), 100 наблюдений')
print('    statistic =', res.statistic)
print('    Есть ли основания отбросить гипотезу?',
      'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')

res = anderson(data2, 'norm')
print('N(0,2), 100 наблюдений')
print('    statistic =', res.statistic)
print('    Есть ли основания отбросить гипотезу?',
      'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')

res = anderson(data3, 'norm')
print('N(0,1), 1000 наблюдений')
print('    statistic =', res.statistic)
print('    Есть ли основания отбросить гипотезу?',
      'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')

res = anderson(data4, 'norm')
print('N(2,1), 1000 наблюдений')
print('    statistic =', res.statistic)
print('    Есть ли основания отбросить гипотезу?',
      'Да' if (res.critical_values[np.where(res.significance_level == 5.)] < res.statistic) else 'Нет')
N(0,1), 100 наблюдений
    statistic = 0.3803493043023565
    Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений
    statistic = 0.5864913565249026
    Есть ли основания отбросить гипотезу? Нет
N(0,1), 1000 наблюдений
    statistic = 0.1886736348968725
    Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений
    statistic = 0.7096445710908483
    Есть ли основания отбросить гипотезу? Нет

R:

In [92]:
%%R

res <- nortest::ad.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- nortest::ad.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- nortest::ad.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- nortest::ad.test(data4)
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------
	Anderson-Darling normality test

data:  data1
A = 0.29819, p-value = 0.5812

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,2), 100 наблюдений ------------------
	Anderson-Darling normality test

data:  data2
A = 0.21992, p-value = 0.8309

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,1), 1000 наблюдений ------------------
	Anderson-Darling normality test

data:  data3
A = 0.25351, p-value = 0.7325

Есть ли основания отбросить гипотезу? Нет 

------------------ N(2,1), 1000 наблюдений ------------------
	Anderson-Darling normality test

data:  data4
A = 0.59499, p-value = 0.1201

Есть ли основания отбросить гипотезу? Нет 

Критерий Крамера - фон Мизеса¶

$$H_0 : F_n(x) = F(x, \theta_0), \quad x \in \mathbb{R},$$

т.е.данные получены из распределения с заданными параметрами.

Положим $F(x, \theta_0) \equiv \Phi(x)$, чтобы проверить, является ли распределение стандартным нормальным.

Python:

In [93]:
res = cramervonmises(data1, 'norm')
print('N(0,1), 100 наблюдений:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')

res = cramervonmises(data2, 'norm')
print('N(0,2), 100 наблюдений:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')

res = cramervonmises(data3, 'norm')
print('N(0,1), 1000 наблюдений:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')

res = cramervonmises(data4, 'norm')
print('N(2,1), 1000 наблюдений:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
    pvalue =  0.36766184664388846
    Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
    pvalue =  0.0007999760848421689
    Есть ли основания отбросить гипотезу? Да 
N(0,1), 1000 наблюдений:
    pvalue =  0.6078851470500846
    Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
    pvalue =  5.3411167155736905e-08
    Есть ли основания отбросить гипотезу? Да 

R:

In [94]:
%%R

res <- goftest::cvm.test(data1, "pnorm")
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- goftest::cvm.test(data2, "pnorm")
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- goftest::cvm.test(data3, "pnorm")
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- goftest::cvm.test(data4, "pnorm")
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	Parameters assumed to be fixed

data:  data1
omega2 = 0.080434, p-value = 0.6904

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,2), 100 наблюдений ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	Parameters assumed to be fixed

data:  data2
omega2 = 1.786, p-value = 3.433e-05

Есть ли основания отбросить гипотезу? Да 

------------------ N(0,1), 1000 наблюдений ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	Parameters assumed to be fixed

data:  data3
omega2 = 0.032382, p-value = 0.9678

Есть ли основания отбросить гипотезу? Нет 

------------------ N(2,1), 1000 наблюдений ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	Parameters assumed to be fixed

data:  data4
omega2 = 226.87, p-value < 2.2e-16

Есть ли основания отбросить гипотезу? Да 

Критерий Колмогорова - Смирнова в модификации Лиллиефорса¶

$$H_0: F_n(x) \in \{ \; F(x,\theta) \; | \; \theta \in \Theta_0 \; \} \equiv $$$$ \equiv \{ \; F(x; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \int^x_{-\infty} e^{-\frac{(t-\mu)^2}{2\sigma^2}} \; | \; \mu \in \mathbb{R}, \sigma > 0 \; \},$$

т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.

Python:

In [95]:
_, pvalue = lilliefors(data1, 'norm')
print('N(0,1), 100 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')

_, pvalue = lilliefors(data2, 'norm')
print('N(0,2), 100 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')

_, pvalue = lilliefors(data3, 'norm')
print('N(0,1), 1000 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')

_, pvalue = lilliefors(data4, 'norm')
print('N(2,1), 1000 наблюдений:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
    pvalue =  0.5508216991827704
    Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
    pvalue =  0.28541751177174657
    Есть ли основания отбросить гипотезу? Нет
N(0,1), 1000 наблюдений:
    pvalue =  0.7354643954471619
    Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
    pvalue =  0.03457540982170773
    Есть ли основания отбросить гипотезу? Да 

R:

In [96]:
%%R

res <- lillie.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- lillie.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- lillie.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- lillie.test(data4)
cat("------------------ N(2,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  data1
D = 0.047378, p-value = 0.8401

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,2), 100 наблюдений ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  data2
D = 0.050275, p-value = 0.7727

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,1), 1000 наблюдений ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  data3
D = 0.017039, p-value = 0.6861

Есть ли основания отбросить гипотезу? Нет 

------------------ N(2,1), 1000 наблюдений ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  data4
D = 0.026296, p-value = 0.09575

Есть ли основания отбросить гипотезу? Нет 

Критерий Шапиро-Франсия¶

$$H_0: F_n(x) \in \{ \; F(x,\theta) \; | \; \theta \in \Theta_0 \; \} \equiv $$$$ \equiv \{ \; F(x; \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \int^x_{-\infty} e^{-\frac{(t-\mu)^2}{2\sigma^2}} \; | \; \mu \in \mathbb{R}, \sigma > 0 \; \},$$

т.е. выборка имеет нормальное распределение, параметры которого заранее не известны.

Критерий Шапиро - Франсия - это упрощение критерия Шапиро - Уилка и для него тоже необходимо, чтобы размер выборки был не более 5000.

Python:

In [97]:
res = shapiroFrancia(data1)
print('N(0,1), 100 наблюдений:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')

res = shapiroFrancia(data2)
print('N(0,2), 100 наблюдений:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')

res = shapiroFrancia(data3)
print('N(0,1), 1000 наблюдений:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')

res = shapiroFrancia(data4)
print('N(2,1), 1000 наблюдений:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < 0.05) else 'Нет')
N(0,1), 100 наблюдений:
    pvalue =  0.5539857967110066
    Есть ли основания отбросить гипотезу? Нет
N(0,2), 100 наблюдений:
    pvalue =  0.030175322690098563
    Есть ли основания отбросить гипотезу? Да 
N(0,1), 1000 наблюдений:
    pvalue =  0.9786190291350362
    Есть ли основания отбросить гипотезу? Нет
N(2,1), 1000 наблюдений:
    pvalue =  0.27727193178255205
    Есть ли основания отбросить гипотезу? Нет

R:

In [98]:
%%R

res <- sf.test(data1)
cat("------------------ N(0,1), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- sf.test(data2)
cat("------------------ N(0,2), 100 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- sf.test(data3)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")

res <- sf.test(data4)
cat("------------------ N(0,1), 1000 наблюдений ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < 0.05) "Да" else "Нет", "\n\n")
------------------ N(0,1), 100 наблюдений ------------------
	Shapiro-Francia normality test

data:  data1
W = 0.98897, p-value = 0.4974

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,2), 100 наблюдений ------------------
	Shapiro-Francia normality test

data:  data2
W = 0.99344, p-value = 0.8462

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,1), 1000 наблюдений ------------------
	Shapiro-Francia normality test

data:  data3
W = 0.99806, p-value = 0.2742

Есть ли основания отбросить гипотезу? Нет 

------------------ N(0,1), 1000 наблюдений ------------------
	Shapiro-Francia normality test

data:  data4
W = 0.99709, p-value = 0.06264

Есть ли основания отбросить гипотезу? Нет 

6. Анализ данных с помощью графиков квантилей, метода огибающих и стандартных процедур проверки гипотез о нормальности (реальные данные)¶

Теперь исследуем теми же методами реальные данные. Рассмотрим следующие выборки:

  • возраст авторов из первого датасета (age, 175 измерений)
    • возраст авторов-женщин (58 измерений)
    • возраст авторов-мужчин (117 измерений)
  • средние оценки пользователей из второго датасета (average_rating, 10000 измерений)
    • оценки для комедий (145 измерений)
    • оценки для драм (406 измерения)

Заметим, что выборки для драм и комедий пересекаются всего по 2 измерениям (см. код ниже).

Визуализация данных¶

Python:

In [99]:
# общие отрезки разбиения для гистограмм age
edges = np.histogram_bin_edges(df1['age'])
In [100]:
plt.figure(figsize=(10, 6))

sns.histplot(data=df1['age'], stat='density', bins=edges, kde=True)

plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('age', fontsize=10, fontweight='bold')
plt.show()
In [101]:
print('Жен:', df1[df1.gender == 'F'].shape[0])
print('Муж:', df1[df1.gender == 'M'].shape[0])
Жен: 58
Муж: 117
In [102]:
plt.figure(figsize=(10, 6))

sns.histplot(data=df1[df1.gender == 'M'].age, stat='density', kde=True,
             bins=edges, alpha=0.3, label='Male', color='blue')
sns.histplot(data=df1[df1.gender == 'F'].age, stat='density', kde=True,
             bins=edges, alpha=0.3, label='Female', color='red')

plt.xticks(fontsize=9)
plt.xlabel('Возраст', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('age (жен и муж)', fontsize=10, fontweight='bold')
plt.legend(title='Пол', title_fontsize=10, fontsize=10)
plt.show()
In [103]:
# общие отрезки разбиения для гистограмм average_rating
# (количество такое, так как average_rating округляется до 2 знаков после запятой)
edges = np.histogram_bin_edges(df2['average_rating'], bins=100)
In [104]:
plt.figure(figsize=(10, 6))

sns.histplot(data=df2['average_rating'], stat='density', bins=edges, kde=True)

plt.xticks(fontsize=9)
plt.xlabel('Средняя оценка', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('average_rating', fontsize=10, fontweight='bold')
plt.show()
In [105]:
df2_explode = df2.explode('genres')

comedy_df = df2_explode[df2_explode.genres == 'Comedy']
print('Комедии:', comedy_df.shape[0])

drama_df = df2_explode[df2_explode.genres == 'Drama']
print('Драмы:', drama_df.shape[0])

dramedy_df = pd.merge(comedy_df, drama_df, how='inner', on=['title'])
print('Драмеди:', dramedy_df.shape[0])
Комедии: 145
Драмы: 406
Драмеди: 2
In [106]:
plt.figure(figsize=(10, 6))

sns.histplot(data=comedy_df['average_rating'], stat='density', kde=True,
             bins=edges, alpha=0.3, label='Comedy')
sns.histplot(data=drama_df['average_rating'], stat='density', kde=True,
             bins=edges, alpha=0.3, label='Drama')

plt.xticks(fontsize=9)
plt.xlabel('Средняя оценка', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('')
plt.title('average_rating (по жанрам)', fontsize=10, fontweight='bold')
plt.legend(title='Жанр', title_fontsize=10, fontsize=10)
plt.show()

R:

In [107]:
%%R

# общие отрезки разбиения для гистограмм age
edges = hist(df1$age, plot = FALSE)
In [108]:
%%R -w 1000 -h 600

hist(df1$age, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Возраст", ylab = "",
     main = "age")

lines(density(df1$age, adjust=), col = "royalblue", lwd = 2)
In [109]:
%%R

cat("Жен:", nrow(df1[df1$gender %like% "F", ]), "\n")
cat("Муж:", nrow(df1[df1$gender %like% "M", ]), "\n")
Жен: 58 
Муж: 117 
In [110]:
%%R -w 1000 -h 600

hist(df1[df1$gender %like% "F", "age"], freq = FALSE,
     breaks = edges$breaks,
     col = my_red,
     xlab = "Возраст", ylab = "",
     main = "age (муж и жен)")
lines(density(df1[df1$gender %like% "F", "age"], adjust=), col = "red", lwd = 3)

hist(df1[df1$gender %like% "M", "age"], freq = FALSE,
     breaks = edges$breaks,
     col = my_blue,
     add = TRUE)
lines(density(df1[df1$gender %like% "M", "age"], adjust=), col = "royalblue", lwd = 3)

legend("topleft", inset=c(0.03, 0), legend=c("Male", "Female"), fill = c(my_blue, my_red), title="Пол")
In [111]:
%%R

# общие отрезки разбиения для гистограмм average_rating
# (количество такое, так как average_rating округляется до 2 знаков после запятой)
edges = hist(df2$average_rating, freq = FALSE, breaks = 100, plot = FALSE)
In [112]:
%%R -w 1000 -h 600

hist(df2$average_rating, freq = FALSE, breaks = edges$breaks,
     col = "lightblue",
     xlab = "Средняя оценка", ylab = "",
     main = "average_rating")

lines(density(df2$average_rating, adjust=), col = "royalblue", lwd = 2)
In [113]:
%%R

comedy_df <- df2[df2$genres %like% "Comedy", ]
cat("Комедии:", nrow(comedy_df), "\n")

drama_df <- df2[df2$genres %like% "Drama", ]
cat("Драмы:", nrow(drama_df), "\n")

dramedy_df <- merge(comedy_df, drama_df, by = "title")
cat("Драмеди:", nrow(dramedy_df), "\n")
Комедии: 145 
Драмы: 406 
Драмеди: 2 
In [114]:
%%R -w 1000 -h 600

my_orange <- t_col("orange", perc = 60)
hist(drama_df$average_rating, freq = FALSE,
     breaks = edges$breaks,
     col = my_orange,
     xlab = "Средняя оценка", ylab = "",
     main = "average_rating (по жанрам)")
lines(density(drama_df$average_rating, adjust=), col = "orange", lwd = 3)

hist(comedy_df$average_rating, freq = FALSE,
     breaks = edges$breaks,
     col = my_blue,
     add = TRUE)
lines(density(comedy_df$average_rating, adjust=), col = "royalblue", lwd = 3)

legend("topleft", inset=c(0.03, 0), legend=c("Comedy", "Drama"), fill = c(my_blue, my_orange), title="Жанр")

Эмпирические функции распределения¶

Будем сравнивать данные с нормальными распределениями, которые имеют соответствующие средние и стандартные отклонения.

Python:

In [115]:
age_mean = statistics.mean(df1.age)
age_std = np.std(df1.age)

female_mean = statistics.mean(df1[df1.gender == 'F'].age)
female_std = np.std(df1[df1.gender == 'F'].age)

male_mean = statistics.mean(df1[df1.gender == 'M'].age)
male_std = np.std(df1[df1.gender == 'M'].age)
In [116]:
fig = plt.figure(figsize=(10,9))
x = np.linspace(0, 100, 1000)

ax1 = plt.subplot(221)
plt.plot(x, stats.norm.cdf(x, age_mean, age_std), color='red')
sns.ecdfplot(data=df1.age)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax1.set_title('age (все)', fontsize=10, fontweight='bold')

ax3 = plt.subplot(223)
plt.plot(x, stats.norm.cdf(x, male_mean, male_std), color='red')
sns.ecdfplot(data=df1[df1.gender == 'F'].age)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax3.set_title('age (жен)', fontsize=10, fontweight='bold')

ax4 = plt.subplot(224)
plt.plot(x, stats.norm.cdf(x, female_mean, female_std), color='red')
sns.ecdfplot(data=df1[df1.gender == 'M'].age)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax4.set_title('age (муж)', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()
In [117]:
all_mean = statistics.mean(df2.average_rating)
all_std = np.std(df2.average_rating)

comedy_mean = statistics.mean(comedy_df.average_rating)
comedy_std = np.std(comedy_df.average_rating)

drama_mean = statistics.mean(drama_df.average_rating)
drama_std = np.std(drama_df.average_rating)
In [118]:
fig = plt.figure(figsize=(10,9))
x = np.linspace(0, 8, 1000)

ax1 = plt.subplot(221)
plt.plot(x, stats.norm.cdf(x, all_mean, all_std), color='red')
sns.ecdfplot(data=df2.average_rating)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax1.set_title('average_rating (все)', fontsize=10, fontweight='bold')

ax3 = plt.subplot(223)
plt.plot(x, stats.norm.cdf(x, comedy_mean, comedy_std), color='red')
sns.ecdfplot(data=comedy_df.average_rating)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax3.set_title('average_rating (комедии)', fontsize=10, fontweight='bold')

ax4 = plt.subplot(224)
plt.plot(x, stats.norm.cdf(x, drama_mean, drama_std), color='red')
sns.ecdfplot(data=drama_df.average_rating)
plt.xticks(fontsize=9)
plt.xlabel('x', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('F(x)', fontsize=10)
ax4.set_title('average_rating (драмы)', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()

R:

In [119]:
%%R

age_mean = mean(df1$age)
age_std = sd(df1$age)

female_mean = mean(df1[df1$gender %like% "F", "age"])
female_std = sd(df1[df1$gender %like% "F", "age"])

male_mean = mean(df1[df1$gender %like% "M", "age"])
male_std = sd(df1[df1$gender %like% "M", "age"])
In [120]:
%%R -w 1000 -h 450

plot(sort(df1$age), pnorm(sort(df1$age), mean = age_mean, sd = age_std), type = "l",
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="age (все)")
plot(ecdf(df1$age), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

par(mfrow = c(1,2))

plot(sort(df1[df1$gender %like% "F", "age"]), pnorm(sort(df1[df1$gender %like% "F", "age"]), mean = female_mean, sd = female_std), type = "l",
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="age (жен)")
plot(ecdf(df1[df1$gender %like% "F", "age"]), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

plot(sort(df1[df1$gender %like% "M", "age"]), pnorm(sort(df1[df1$gender %like% "M", "age"]), mean = male_mean, sd = male_std), type = "l",
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="age (муж)")
plot(ecdf(df1[df1$gender %like% "M", "age"]), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)
In [121]:
%%R

all_mean = mean(df2$average_rating)
all_std = sd(df2$average_rating)

comedy_mean = mean(comedy_df$average_rating)
comedy_std = sd(comedy_df$average_rating)

drama_mean = mean(drama_df$average_rating)
drama_std = sd(drama_df$average_rating)
In [122]:
%%R -w 1000 -h 450

plot(sort(df2$average_rating), pnorm(sort(df2$average_rating), mean = all_mean, sd = all_std), type = "l",
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="average_rating (все)")
plot(ecdf(df2$average_rating), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

par(mfrow = c(1,2))

plot(sort(comedy_df$average_rating), pnorm(sort(comedy_df$average_rating), mean = comedy_mean, sd = comedy_std), type = "l",
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="average_rating (комедии)")
plot(ecdf(comedy_df$average_rating), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

plot(sort(drama_df$average_rating), pnorm(sort(drama_df$average_rating), mean = drama_mean, sd = drama_std), type = "l",
     lwd = 2, col = "red",
     xlab = "x", ylab = "F(x)",
     main="average_rating (драмы)")
plot(ecdf(drama_df$average_rating), pch=20, cex=0, col="royalblue", lwd=4, add=TRUE)

Квантили¶

Python:

In [123]:
fig = plt.figure(figsize=(10,9))

ax1 = plt.subplot(221)
sm.qqplot(df1.age, ax=ax1, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax1.set_title('age (все)', fontsize=10, fontweight='bold')

ax3 = plt.subplot(223)
sm.qqplot(df1[df1.gender == 'F'].age, ax=ax3, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax3.set_title('age (жен)', fontsize=10, fontweight='bold')

ax4 = plt.subplot(224)
sm.qqplot(df1[df1.gender == 'M'].age, ax=ax4, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax4.set_title('age (муж)', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()
In [124]:
fig = plt.figure(figsize=(10,9))

ax1 = plt.subplot(221)
sm.qqplot(df2.average_rating, ax=ax1, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax1.set_title('average_rating (все)', fontsize=10, fontweight='bold')

ax3 = plt.subplot(223)
sm.qqplot(comedy_df.average_rating, ax=ax3, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax3.set_title('average_rating (комедии)', fontsize=10, fontweight='bold')

ax4 = plt.subplot(224)
sm.qqplot(drama_df.average_rating, ax=ax4, line='s')
plt.xticks(fontsize=9)
plt.xlabel('Теоретические квантили', fontsize=10)
plt.yticks(fontsize=9)
plt.ylabel('Выборочные квантили', fontsize=10)
ax4.set_title('average_rating (драмы)', fontsize=10, fontweight='bold')

fig.tight_layout(h_pad=2, w_pad=2)

plt.show()

R:

In [125]:
%%R -w 1000 -h 450

qqnorm(df1$age,
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "age (все)")
qqline(df1$age,
       col = "red", lwd = 2)

par(mfrow = c(1,2))

qqnorm(df1[df1$gender %like% "F", "age"],
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "age (жен)")
qqline(df1[df1$gender %like% "F", "age"],
       col = "red", lwd = 2)

qqnorm(df1[df1$gender %like% "M", "age"],
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "age (муж)")
qqline(df1[df1$gender %like% "M", "age"],
       col = "red", lwd = 2)
In [126]:
%%R -w 1000 -h 450

qqnorm(df2$average_rating,
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "average_rating (все)")
qqline(df2$average_rating,
       col = "red", lwd = 2)

par(mfrow = c(1,2))

qqnorm(comedy_df$average_rating,
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "average_rating (комедии)")
qqline(comedy_df$average_rating,
       col = "red", lwd = 2)

qqnorm(drama_df$average_rating,
       col = "royalblue", pch = 19,
       xlab = "Теоретические квантили", ylab = "Выборочные квантили",
       main = "average_rating (драмы)")
qqline(drama_df$average_rating,
       col = "red", lwd = 2)

Метод огибающих¶

Python:

In [127]:
fig, ax = plt.subplots(nrows=2, ncols=2, figsize = (10, 9))

envelmet(df1.age, 'age (все)', ax[0, 0])
fig.delaxes(ax[0, 1])
envelmet(df1[df1.gender == 'F'].age, 'age (жен)', ax[1, 0])
envelmet(df1[df1.gender == 'M'].age, 'age (муж)', ax[1, 1])

fig.tight_layout()

plt.show()
In [128]:
fig, ax = plt.subplots(nrows=2, ncols=2, figsize = (10, 9))

envelmet(df2.average_rating, 'average_rating (все)', ax[0, 0])
fig.delaxes(ax[0, 1])
envelmet(comedy_df.average_rating, 'average_rating (комедии)', ax[1, 0])
envelmet(drama_df.average_rating, 'average_rating (драмы)', ax[1, 1])

fig.tight_layout()

plt.show()

R:

In [129]:
%%R -w 1000 -h 450

envelmet(df1$age, "age (все)")

par(mfrow = c(1,2))

envelmet(df1[df1$gender %like% "F", "age"], "age (жен)")
envelmet(df1[df1$gender %like% "M", "age"], "age (муж)")
In [130]:
%%R -w 1000 -h 450

envelmet(df2$average_rating, "average_rating (все)")

par(mfrow = c(1,2))

envelmet(comedy_df$average_rating, "average_rating (комедии)")
envelmet(drama_df$average_rating, "average_rating (драмы)")

Стандартные процедуры проверки нормальности¶

Критерий Колмогорова - Смирнова¶

Воспользуемся двухвыборочным критерием.

Есть 2 выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m$ с функциями распределения $F(x)$ и $G(x)$ соответственно. $$H_0: F(x) \equiv G(x),$$ т.е. выборки принадлежат к одному распределению.

Если полученный при использовании критерия $pvalue$ меньше заданного заранее уровня значимости $\alpha$ (положим $\alpha = 0.5$), то $H_0$ должна быть отвергнута, иначе - нет оснований отбросить нулевую гипотезу.

Во-первых, чтобы проверить нормальность данных, семплируем по выборке из нормальных распределений со средними и стандартными отклонениями, соответствующими этим же величинам у исследуемых распределений.

Во-вторых, проверим, одинаково ли распределены данные:

  • по мужчинам и женщинам
  • по комедиям и драмам

Python:

In [131]:
alpha = 0.05
In [132]:
age_norm = np.random.normal(age_mean, age_std, size=1000)
female_norm = np.random.normal(female_mean, female_std, size=100)
male_norm = np.random.normal(male_mean, male_mean, size=100)
In [133]:
_, pvalue = kstest(df1.age, age_norm)
print('age - нормальное распределение:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')

_, pvalue = kstest(df1[df1.gender == 'F'].age, female_norm)
print('age (жен) - нормальное распределение:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')

_, pvalue = kstest(df1[df1.gender == 'M'].age, male_norm)
print('age (муж) - нормальное распределение:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
age - нормальное распределение:
    pvalue = 0.12182158757492134
    Есть ли основания отбросить гипотезу? Нет
age (жен) - нормальное распределение:
    pvalue = 0.910886686614063
    Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
    pvalue = 6.907798016052148e-12
    Есть ли основания отбросить гипотезу? Да
In [134]:
_, pvalue = kstest(df1[df1.gender == 'F'].age, df1[df1.gender == 'M'].age)
print('Данные по женщинам и мужчинам одинаково распределены:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
Данные по женщинам и мужчинам одинаково распределены:
    pvalue = 0.7227115888214679
    Есть ли основания отбросить гипотезу? Нет
In [135]:
all_norm = np.random.normal(all_mean, all_std, size=1000)
comedy_norm = np.random.normal(comedy_mean, comedy_std, size=100)
drama_norm = np.random.normal(drama_mean, comedy_mean, size=100)
In [136]:
_, pvalue = kstest(df2.average_rating, all_norm)
print('average_rating - нормальное распределение:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')

_, pvalue = kstest(comedy_df.average_rating, comedy_norm)
print('average_rating (комедии) - нормальное распределение:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')

_, pvalue = kstest(drama_df.average_rating, drama_norm)
print('average_rating (драмы) - нормальное распределение:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
    pvalue = 0.001585106579657839
    Есть ли основания отбросить гипотезу? Да
average_rating (комедии) - нормальное распределение:
    pvalue = 0.8805561402063771
    Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
    pvalue = 8.713468542532803e-14
    Есть ли основания отбросить гипотезу? Да
In [137]:
_, pvalue = kstest(comedy_df.average_rating, drama_df.average_rating)
print('Данные по комедиям и драмам одинаково распределены:')
print('    pvalue =', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да' if (pvalue < alpha) else 'Нет')
Данные по комедиям и драмам одинаково распределены:
    pvalue = 0.020898085404186383
    Есть ли основания отбросить гипотезу? Да

R:

In [138]:
%%R

alpha <- 0.05
In [139]:
%%R

age_norm <- rnorm(1000, mean = age_mean, sd = age_std)
female_norm <- rnorm(100, mean = female_mean, sd = female_std)
male_norm <- rnorm(100, mean = male_mean, sd = male_mean)
In [140]:
%%R

res <- ks.test(df1$age, age_norm)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- ks.test(df1[df1$gender %like% "F", "age"], female_norm)
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- ks.test(df1[df1$gender %like% "M", "age"], male_norm)
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------
	Asymptotic two-sample Kolmogorov-Smirnov test

data:  df1$age and age_norm
D = 0.11271, p-value = 0.04545
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Да 

------------------ age (жен) - нормальное распределение ------------------
	Exact two-sample Kolmogorov-Smirnov test

data:  df1[df1$gender %like% "F", "age"] and female_norm
D = 0.20379, p-value = 0.07555
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Нет 

------------------ age (муж) - нормальное распределение ------------------
	Asymptotic two-sample Kolmogorov-Smirnov test

data:  df1[df1$gender %like% "M", "age"] and male_norm
D = 0.42291, p-value = 8.418e-09
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Да 

In [141]:
%%R

res <- ks.test(df1[df1$gender %like% "F", "age"], df1[df1$gender %like% "M", "age"])
cat("------------------ Данные по женщинам и мужчинам одинаково распределены ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ Данные по женщинам и мужчинам одинаково распределены ------------------
	Exact two-sample Kolmogorov-Smirnov test

data:  df1[df1$gender %like% "F", "age"] and df1[df1$gender %like% "M", "age"]
D = 0.1064, p-value = 0.6048
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Нет 

In [142]:
%%R

all_norm <- rnorm(1000, mean = all_mean, sd = all_std)
comedy_norm <- rnorm(100, mean = comedy_mean, sd = comedy_std)
drama_norm <- rnorm(100, mean = drama_mean, sd = comedy_mean)
In [143]:
%%R

res <- ks.test(df2$average_rating, all_norm)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- ks.test(comedy_df$average_rating, comedy_norm)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- ks.test(drama_df$average_rating, drama_norm)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------
	Asymptotic two-sample Kolmogorov-Smirnov test

data:  df2$average_rating and all_norm
D = 0.0667, p-value = 0.0006139
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Да 

------------------ average_rating (комедии) - нормальное распределение ------------------
	Asymptotic two-sample Kolmogorov-Smirnov test

data:  comedy_df$average_rating and comedy_norm
D = 0.10724, p-value = 0.504
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Нет 

------------------ average_rating (драмы) - нормальное распределение ------------------
	Asymptotic two-sample Kolmogorov-Smirnov test

data:  drama_df$average_rating and drama_norm
D = 0.43015, p-value = 2.547e-13
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Да 

In [144]:
%%R

res <- ks.test(comedy_df$average_rating, drama_df$average_rating)
cat("------------------ Данные по комедиям и драмам одинаково распределены ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ Данные по комедиям и драмам одинаково распределены ------------------
	Asymptotic two-sample Kolmogorov-Smirnov test

data:  comedy_df$average_rating and drama_df$average_rating
D = 0.14433, p-value = 0.02332
alternative hypothesis: two-sided

Есть ли основания отбросить гипотезу? Да 

Критерий Шапиро - Уилка¶

Так как для критерия Шапиро-Уилка необходимо, чтобы размер выборки был не более 5000, возьмём только половину от всех данных по рейтингам (отсортируем по возрастанию average_rating и возьмём каждое второе значение).

Python:

In [145]:
_, pvalue = shapiro(df1.age)
print('age - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = shapiro(df1[df1.gender == 'F'].age)
print('age (жен) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = shapiro(df1[df1.gender == 'M'].age)
print('age (муж) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
age - нормальное распределение:
    pvalue =  0.0003017295675817877
    Есть ли основания отбросить гипотезу? Да 
age (жен) - нормальное распределение:
    pvalue =  0.24597471952438354
    Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
    pvalue =  0.0003075910790357739
    Есть ли основания отбросить гипотезу? Да 
In [146]:
df2_half = df2.sort_values('average_rating').reset_index().iloc[lambda x: x.index % 2 == 0].reset_index()
In [147]:
_, pvalue = shapiro(df2_half.average_rating)
print('average_rating - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = shapiro(comedy_df.average_rating)
print('average_rating (комедии) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = shapiro(drama_df.average_rating)
print('average_rating (драмы) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
    pvalue =  0.0
    Есть ли основания отбросить гипотезу? Да 
average_rating (комедии) - нормальное распределение:
    pvalue =  0.7259739637374878
    Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
    pvalue =  3.4041782726035308e-09
    Есть ли основания отбросить гипотезу? Да 

R:

In [148]:
%%R

res <- shapiro.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- shapiro.test(df1[df1$gender %like% "F", "age"])
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- shapiro.test(df1[df1$gender %like% "M", "age"])
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------
	Shapiro-Wilk normality test

data:  df1$age
W = 0.96624, p-value = 0.0003017

Есть ли основания отбросить гипотезу? Да 

------------------ age (жен) - нормальное распределение ------------------
	Shapiro-Wilk normality test

data:  df1[df1$gender %like% "F", "age"]
W = 0.97398, p-value = 0.246

Есть ли основания отбросить гипотезу? Нет 

------------------ age (муж) - нормальное распределение ------------------
	Shapiro-Wilk normality test

data:  df1[df1$gender %like% "M", "age"]
W = 0.95093, p-value = 0.0003076

Есть ли основания отбросить гипотезу? Да 

In [149]:
%%R

df2_half <- df2 %>%
  arrange(average_rating) %>%
  mutate(index = row_number()) %>%
  filter(index %% 2 == 0) %>%
  select(-index) %>%
  mutate(index = row_number()) %>%
  select(-index)
In [150]:
%%R

res <- shapiro.test(df2_half$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- shapiro.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- shapiro.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------
	Shapiro-Wilk normality test

data:  df2_half$average_rating
W = 0.91881, p-value < 2.2e-16

Есть ли основания отбросить гипотезу? Да 

------------------ average_rating (комедии) - нормальное распределение ------------------
	Shapiro-Wilk normality test

data:  comedy_df$average_rating
W = 0.99322, p-value = 0.726

Есть ли основания отбросить гипотезу? Нет 

------------------ average_rating (драмы) - нормальное распределение ------------------
	Shapiro-Wilk normality test

data:  drama_df$average_rating
W = 0.95912, p-value = 3.405e-09

Есть ли основания отбросить гипотезу? Да 

Критерий Андерсона - Дарлинга¶

Python:


На Python есть реализация критерия Андерсона для $k$ выборок. С его помощью (также как с двухвыборочным критерием Колмогорова - Смирнова) проверим и то, что данные нормальны, и то, что данные по комедиям и драмам распределены одинаковы.

In [151]:
_, _, pvalue = anderson_ksamp([df1.age, age_norm])
print('age - нормальное распределеныие:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, _, pvalue = anderson_ksamp([df1[df1.gender == 'F'].age, female_norm])
print('age (жен) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, _, pvalue = anderson_ksamp([df1[df1.gender == 'M'].age, male_norm])
print('age (муж) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
age - нормальное распределеныие:
    pvalue =  0.25
    Есть ли основания отбросить гипотезу? Нет
age (жен) - нормальное распределение:
    pvalue =  0.25
    Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
    pvalue =  0.001
    Есть ли основания отбросить гипотезу? Да 
In [152]:
_, _, pvalue = anderson_ksamp([df1[df1.gender == 'F'].age, df1[df1.gender == 'M'].age])
print('Данные по женщинам и мужчинам одинаково распределены:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
Данные по женщинам и мужчинам одинаково распределены:
    pvalue =  0.25
    Есть ли основания отбросить гипотезу? Нет
In [153]:
_, _, pvalue = anderson_ksamp([df2.average_rating, all_norm])
print('average_rating - нормальное распределеныие:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, _, pvalue = anderson_ksamp([comedy_df.average_rating, comedy_norm])
print('average_rating (комедии) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, _, pvalue = anderson_ksamp([drama_df.average_rating, drama_norm])
print('average_rating (драмы) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределеныие:
    pvalue =  0.0033431320434570704
    Есть ли основания отбросить гипотезу? Да 
average_rating (комедии) - нормальное распределение:
    pvalue =  0.25
    Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
    pvalue =  0.001
    Есть ли основания отбросить гипотезу? Да 
In [154]:
_, _, pvalue = anderson_ksamp([comedy_df.average_rating, drama_df.average_rating])
print('Данные по комедиям и драмам одинаково распределены:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
Данные по комедиям и драмам одинаково распределены:
    pvalue =  0.03658924476599451
    Есть ли основания отбросить гипотезу? Да 

R:

In [155]:
%%R

res <- nortest::ad.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- nortest::ad.test(df1[df1$gender %like% "F", "age"])
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- nortest::ad.test(df1[df1$gender %like% "M", "age"])
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------
	Anderson-Darling normality test

data:  df1$age
A = 1.3089, p-value = 0.002072

Есть ли основания отбросить гипотезу? Да 

------------------ age (жен) - нормальное распределение ------------------
	Anderson-Darling normality test

data:  df1[df1$gender %like% "F", "age"]
A = 0.5091, p-value = 0.1906

Есть ли основания отбросить гипотезу? Нет 

------------------ age (муж) - нормальное распределение ------------------
	Anderson-Darling normality test

data:  df1[df1$gender %like% "M", "age"]
A = 1.0999, p-value = 0.006709

Есть ли основания отбросить гипотезу? Да 

In [156]:
%%R

res <- nortest::ad.test(df2$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- nortest::ad.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- nortest::ad.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------
	Anderson-Darling normality test

data:  df2$average_rating
A = 57.075, p-value < 2.2e-16

Есть ли основания отбросить гипотезу? Да 

------------------ average_rating (комедии) - нормальное распределение ------------------
	Anderson-Darling normality test

data:  comedy_df$average_rating
A = 0.30901, p-value = 0.5541

Есть ли основания отбросить гипотезу? Нет 

------------------ average_rating (драмы) - нормальное распределение ------------------
	Anderson-Darling normality test

data:  drama_df$average_rating
A = 4.5655, p-value = 2.455e-11

Есть ли основания отбросить гипотезу? Да 

Критерий Крамера - фон Мизеса¶

Python:


На Python есть реализация критерия Крамера - фон Мизеса для $k$ выборок. С его помощью (также как и с двухвыборочным критерием Колмогорова - Смирнова) проверим и то, что данные нормальны, и то, что данные по мужчинам/женщинам и комедиям/драмам распределены одинаковы.

In [157]:
res = cramervonmises_2samp(df1.age, age_norm)
print('age - нормальное распределение:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')

res = cramervonmises_2samp(df1[df1.gender == 'F'].age, female_norm)
print('age (жен) - нормальное распределение:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')

res = cramervonmises_2samp(df1[df1.gender == 'M'].age, male_norm)
print('age (муж) - нормальное распределение:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
age - нормальное распределение:
    pvalue =  0.39133160782598164
    Есть ли основания отбросить гипотезу? Нет
age (жен) - нормальное распределение:
    pvalue =  0.9345779001744031
    Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
    pvalue =  1.1599448268562185e-07
    Есть ли основания отбросить гипотезу? Да 
In [158]:
res = cramervonmises_2samp(df1[df1.gender == 'F'].age, df1[df1.gender == 'M'].age)
print('Данные по женщинам и мужчинам одинаково распределены:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
Данные по женщинам и мужчинам одинаково распределены:
    pvalue =  0.5198899109819444
    Есть ли основания отбросить гипотезу? Нет
In [159]:
res = cramervonmises_2samp(df2.average_rating, all_norm)
print('average_rating - нормальное распределение:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')

res = cramervonmises_2samp(comedy_df.average_rating, comedy_norm)
print('average_rating (комедии) - нормальное распределение:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')

res = cramervonmises_2samp(drama_df.average_rating, drama_norm)
print('average_rating (драмы) - нормальное распределение:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
    pvalue =  0.00012693698164067957
    Есть ли основания отбросить гипотезу? Да 
average_rating (комедии) - нормальное распределение:
    pvalue =  0.8909349351026549
    Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
    pvalue =  2.1094737068239056e-10
    Есть ли основания отбросить гипотезу? Да 
In [160]:
res = cramervonmises_2samp(comedy_df.average_rating, drama_df.average_rating)
print('Данные по комедиям и драмам одинаково распределены:')
print('    pvalue = ', res.pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res.pvalue < alpha) else 'Нет')
Данные по комедиям и драмам одинаково распределены:
    pvalue =  0.03926253837913829
    Есть ли основания отбросить гипотезу? Да 

R:

In [161]:
%%R

res <- goftest::cvm.test(df1$age, "pnorm", mean = age_mean, sd = age_std)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- goftest::cvm.test(df1[df1$gender %like% 'F', "age"], "pnorm", mean = female_mean, sd = female_std)
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- goftest::cvm.test(df1[df1$gender %like% 'M', "age"], "pnorm", mean = male_mean, sd = male_std)
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	with parameters mean = 42.8457142857143, sd = 10.9052434624528
	Parameters assumed to be fixed

data:  df1$age
omega2 = 0.19279, p-value = 0.282

Есть ли основания отбросить гипотезу? Нет 

------------------ age (жен) - нормальное распределение ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	with parameters mean = 41.6206896551724, sd = 10.8786408887419
	Parameters assumed to be fixed

data:  df1[df1$gender %like% "F", "age"]
omega2 = 0.084203, p-value = 0.6697

Есть ли основания отбросить гипотезу? Нет 

------------------ age (муж) - нормальное распределение ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	with parameters mean = 43.4529914529915, sd = 10.913844202865
	Parameters assumed to be fixed

data:  df1[df1$gender %like% "M", "age"]
omega2 = 0.14668, p-value = 0.4003

Есть ли основания отбросить гипотезу? Нет 

In [162]:
%%R

res <- goftest::cvm.test(df2$average_rating, "pnorm", mean = all_mean, sd = all_std)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- goftest::cvm.test(comedy_df$average_rating, "pnorm", mean = comedy_mean, sd = comedy_std)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- goftest::cvm.test(drama_df$average_rating, "pnorm", mean = drama_mean, sd = drama_std)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	with parameters mean = 4.068577, sd = 0.335358601515096
	Parameters assumed to be fixed

data:  df2$average_rating
omega2 = 8.5192, p-value < 2.2e-16

Есть ли основания отбросить гипотезу? Да 

------------------ average_rating (комедии) - нормальное распределение ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	with parameters mean = 4.02427586206897, sd = 0.266092467291777
	Parameters assumed to be fixed

data:  comedy_df$average_rating
omega2 = 0.038392, p-value = 0.9418

Есть ли основания отбросить гипотезу? Нет 

------------------ average_rating (драмы) - нормальное распределение ------------------
	Cramer-von Mises test of goodness-of-fit
	Null hypothesis: Normal distribution
	with parameters mean = 4.02997536945813, sd = 0.294512360107409
	Parameters assumed to be fixed

data:  drama_df$average_rating
omega2 = 0.80123, p-value = 0.007232

Есть ли основания отбросить гипотезу? Да 

Критерий Колмогора - Смирнова в модификации Лиллиефорса¶

Python:

In [163]:
_, pvalue = lilliefors(df1.age, 'norm')
print('age - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = lilliefors(df1[df1.gender == 'F'].age, 'norm')
print('age (жен) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = lilliefors(df1[df1.gender == 'M'].age, 'norm')
print('age (муж) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
age - нормальное распределение:
    pvalue =  0.001761024787129434
    Есть ли основания отбросить гипотезу? Да 
age (жен) - нормальное распределение:
    pvalue =  0.26096030891637234
    Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
    pvalue =  0.0206531490309459
    Есть ли основания отбросить гипотезу? Да 
In [164]:
_, pvalue = lilliefors(df2.average_rating, 'norm')
print('average_rating - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = lilliefors(comedy_df.average_rating, 'norm')
print('average_rating (комедии) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')

_, pvalue = lilliefors(drama_df.average_rating, 'norm')
print('average_rating (драмы) - нормальное распределение:')
print('    pvalue = ', pvalue)
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (pvalue < alpha) else 'Нет')
average_rating - нормальное распределение:
    pvalue =  0.0009999999999998899
    Есть ли основания отбросить гипотезу? Да 
average_rating (комедии) - нормальное распределение:
    pvalue =  0.9073222722620489
    Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
    pvalue =  0.0009999999999998899
    Есть ли основания отбросить гипотезу? Да 

R:

In [165]:
%%R

res <- lillie.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- lillie.test(df1[df1$gender %like% 'F', "age"])
cat("------------------ age (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- lillie.test(df1[df1$gender %like% 'M', "age"])
cat("------------------ age (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  df1$age
D = 0.092914, p-value = 0.0008559

Есть ли основания отбросить гипотезу? Да 

------------------ age (комедии) - нормальное распределение ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  df1[df1$gender %like% "F", "age"]
D = 0.093698, p-value = 0.2327

Есть ли основания отбросить гипотезу? Нет 

------------------ age (драмы) - нормальное распределение ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  df1[df1$gender %like% "M", "age"]
D = 0.09319, p-value = 0.01424

Есть ли основания отбросить гипотезу? Да 

In [166]:
%%R

res <- lillie.test(df2$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- lillie.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- lillie.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  df2$average_rating
D = 0.053433, p-value < 2.2e-16

Есть ли основания отбросить гипотезу? Да 

------------------ average_rating (комедии) - нормальное распределение ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  comedy_df$average_rating
D = 0.038619, p-value = 0.8609

Есть ли основания отбросить гипотезу? Нет 

------------------ average_rating (драмы) - нормальное распределение ------------------
	Lilliefors (Kolmogorov-Smirnov) normality test

data:  drama_df$average_rating
D = 0.11116, p-value = 4.371e-13

Есть ли основания отбросить гипотезу? Да 

Критерий Шапиро - Франсия¶

Критерий Шапиро - Франсия - это упрощение критерия Шапиро - Уилка и для него тоже необходимо, чтобы размер выборки был не более 5000. Опять берём только половину данных по рейтингам.

Python:

In [167]:
res = shapiroFrancia(df1.first_published)
print('age - нормальное распределение:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')

res = shapiroFrancia(df1[df1.gender == 'F'].age)
print('age (жен) - нормальное распределение:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')

res = shapiroFrancia(df1[df1.gender == 'M'].age)
print('age (муж) - нормальное распределение:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
age - нормальное распределение:
    pvalue =  1.750186564963448e-18
    Есть ли основания отбросить гипотезу? Да 
age (жен) - нормальное распределение:
    pvalue =  0.1360614634580179
    Есть ли основания отбросить гипотезу? Нет
age (муж) - нормальное распределение:
    pvalue =  0.0005470998390365658
    Есть ли основания отбросить гипотезу? Да 
In [168]:
res = shapiroFrancia(df2_half.average_rating)
print('average_rating - нормальное распределение:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')

res = shapiroFrancia(comedy_df.average_rating)
print('average_rating (комедии) - нормальное распределение:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')

res = shapiroFrancia(drama_df.average_rating)
print('average_rating (драмы) - нормальное распределение:')
print('    pvalue = ', res['p-value'])
print('    Есть ли основания отбросить гипотезу?', 'Да ' if (res['p-value'] < alpha) else 'Нет')
average_rating - нормальное распределение:
    pvalue =  3.293768000557143e-46
    Есть ли основания отбросить гипотезу? Да 
average_rating (комедии) - нормальное распределение:
    pvalue =  0.46524972040289425
    Есть ли основания отбросить гипотезу? Нет
average_rating (драмы) - нормальное распределение:
    pvalue =  2.0959308700711756e-08
    Есть ли основания отбросить гипотезу? Да 

R:

In [169]:
%%R

res <- sf.test(df1$age)
cat("------------------ age - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- sf.test(df1[df1$gender %like% "F", "age"])
cat("------------------ age (жен) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- sf.test(df1[df1$gender %like% "M", "age"])
cat("------------------ age (муж) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ age - нормальное распределение ------------------
	Shapiro-Francia normality test

data:  df1$age
W = 0.96513, p-value = 0.0004398

Есть ли основания отбросить гипотезу? Да 

------------------ age (жен) - нормальное распределение ------------------
	Shapiro-Francia normality test

data:  df1[df1$gender %like% "F", "age"]
W = 0.96964, p-value = 0.1361

Есть ли основания отбросить гипотезу? Нет 

------------------ age (муж) - нормальное распределение ------------------
	Shapiro-Francia normality test

data:  df1[df1$gender %like% "M", "age"]
W = 0.95049, p-value = 0.0005471

Есть ли основания отбросить гипотезу? Да 

In [170]:
%%R

res <- sf.test(df2_half$average_rating)
cat("------------------ average_rating - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- sf.test(comedy_df$average_rating)
cat("------------------ average_rating (комедии) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")

res <- sf.test(drama_df$average_rating)
cat("------------------ average_rating (драмы) - нормальное распределение ------------------")
print(res)
cat("Есть ли основания отбросить гипотезу?",
    if (res$p.value < alpha) "Да" else "Нет", "\n\n")
------------------ average_rating - нормальное распределение ------------------
	Shapiro-Francia normality test

data:  df2_half$average_rating
W = 0.91803, p-value < 2.2e-16

Есть ли основания отбросить гипотезу? Да 

------------------ average_rating (комедии) - нормальное распределение ------------------
	Shapiro-Francia normality test

data:  comedy_df$average_rating
W = 0.99158, p-value = 0.4652

Есть ли основания отбросить гипотезу? Нет 

------------------ average_rating (драмы) - нормальное распределение ------------------
	Shapiro-Francia normality test

data:  drama_df$average_rating
W = 0.95811, p-value = 2.096e-08

Есть ли основания отбросить гипотезу? Да 

Заметим:

  1. В первом датасете:

    • Данные по всем возрастам прошли половину критериев на нормальность
    • Данные по женщинам можем считать нормальными (прошли все критерии)
    • Данные по мужчинам прошли только по критерию Крамера фон Мизеса
    • Данные по женщинам и мужчинам распределены одинаково (прошли все критерии)

    При доверительном уровне $\alpha = 0.1$ картина не меняется. А при доверительном уровне $\alpha = 0.01$ данные по мужчинам проходят ещё и критерий Лиллиефорса.

    Поэтому будем считать, что распределения данных по всем возрастам и по возрастам мужчин достаточно близки к нормальным, чтобы можно было пренебречь этой разницей (скорее всего, у этих распределений более тяжёлые хвосты: это видно из гистограммы; к тому же известно, что сильнее критерия Колмогорова - Смирнова критерий Крамера - фон Мизеса реагирует на отклонения в среднем, а критерий Андерсона - Дарлинга - на тяжёлые хвосты).

  2. Во втором датасете:

    • Нормально распределёнными можно считать только данные по комедиям
    • Данные по комедиям и драмам не распределены одинаково

7. Применение критериев Стьюдента, Уилкоксона - Манна - Уитни, Фишера, Левене, Бартлетта и Флигнера - Килина¶

Критерии Стьюдента¶

Для каждого из критериев формулируется нулевая гипотеза $H_0$ и альтернативная гипотеза. Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается в пользу альтернативы, иначе - нет оснований отбросить гипотезу.

Также все разновидности критерия Стьюдента основаны на дополнительном предположении о нормальности выборки данных, поэтому будем использовать эти разновидности на данных о возрастах авторов, которые считаем нормально распределёнными.

Одновыборочный критерий Стьюдента¶

Дана выборка $X_1, \dots, X_n \sim N(\mu, \alpha).$ $$H_0: \mathbb{E}X = \mu_0,$$ т.е. матожидание выборки равно заданному числу.

Альтернативные гипотезы: $$H_1: \mathbb{E}X \neq \mu_0 \quad (\text{двухсторонняя гипотеза})$$ $$H_1': \mathbb{E}X > \mu_0 \quad (\text{односторонняя гипотеза})$$ $$H_1'': \mathbb{E}X < \mu_0 \quad (\text{односторонняя гипотеза})$$

Проверим, равен ли средний возраст авторов 40 годам.

Python:

In [171]:
print('mean(X) =', statistics.mean(df1.age), 'vs mu0 = 40')

print('E(X) != mu0')
res = ttest(df1.age, 40, alternative = 'two-sided', confidence=0.9)
print('    0.90:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'two-sided', confidence=0.95)
print('    0.95:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'two-sided', confidence=0.99)
print('    0.99:   pvalue =', res['p-val'][0], '  power =', res['power'][0])

print('E(X) > mu0')
res = ttest(df1.age, 40, alternative = 'greater', confidence=0.9)
print('    0.90:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'greater', confidence=0.95)
print('    0.95:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'greater', confidence=0.99)
print('    0.99:   pvalue =', res['p-val'][0], '  power =', res['power'][0])

print('E(X) < mu0')
res = ttest(df1.age, 40, alternative = 'less', confidence=0.9)
print('    0.90:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'less', confidence=0.95)
print('    0.95:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1.age, 40, alternative = 'less', confidence=0.99)
print('    0.99:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
mean(X) = 42.84571428571429 vs mu0 = 40
E(X) != mu0
    0.90:   pvalue = 0.0006984719380686418   power = 0.9296182976236684
    0.95:   pvalue = 0.0006984719380686418   power = 0.9296182976236684
    0.99:   pvalue = 0.0006984719380686418   power = 0.9296182976236684
E(X) > mu0
    0.90:   pvalue = 0.0003492359690343209   power = 0.9635704269241538
    0.95:   pvalue = 0.0003492359690343209   power = 0.9635704269241538
    0.99:   pvalue = 0.0003492359690343209   power = 0.9635704269241538
E(X) < mu0
    0.90:   pvalue = 0.9996507640309656   power = 1.8520793543252978e-07
    0.95:   pvalue = 0.9996507640309656   power = 1.8520793543252978e-07
    0.99:   pvalue = 0.9996507640309656   power = 1.8520793543252978e-07

R:

In [172]:
%%R

cat("mean(X) =", mean(df1$age), "vs mu0 = 40")

cat("\nE(X) == mu0\n")
res = t.test(df1$age, mu = 40, alternative = "two.sided", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "two.sided", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "two.sided", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, "\n")

cat("\nE(X) < mu0\n")
res = t.test(df1$age, mu = 40, alternative = "greater", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "greater", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "greater", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, "\n")

cat("\nE(X) > mu0\n")
res = t.test(df1$age, mu = 40, alternative = "less", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "less", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, "\n")
res = t.test(df1$age, mu = 40, alternative = "less", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, "\n")
mean(X) = 42.84571 vs mu0 = 40
E(X) == mu0
    0.90:  pvalue = 0.0006984719 
    0.95:  pvalue = 0.0006984719 
    0.99:  pvalue = 0.0006984719 

E(X) < mu0
    0.90:  pvalue = 0.000349236 
    0.95:  pvalue = 0.000349236 
    0.99:  pvalue = 0.000349236 

E(X) > mu0
    0.90:  pvalue = 0.9996508 
    0.95:  pvalue = 0.9996508 
    0.99:  pvalue = 0.9996508 

Средний возраст авторов больше 40.

Двухвыборочный критерий Стьюдента¶

Даны независимые выборки $X_1, \dots, X_n \sim N(\mu_1, \alpha_1)$ и $Y_1, \dots, Y_m \sim N(\mu_2, \alpha_2).$ $$H_0: \mathbb{E}X = \mathbb{E}Y,$$ т.е. матожидания выборок равны.

Альтернативные гипотезы: $$H_1: \mathbb{E}X \neq \mathbb{E}Y \quad (\text{двухсторонняя гипотеза})$$ $$H_1': \mathbb{E}X > \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$ $$H_1'': \mathbb{E}X < \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$

Проверим, равны ли средний возраст авторов-мужчин и авторов-женщин.

Python:

In [173]:
print('mean(X) =', statistics.mean(df1[df1.gender == 'M'].age),
      'vs mean(Y) =', statistics.mean(df1[df1.gender == 'F'].age))

print('E(X) != E(Y)')
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'two-sided', confidence=0.9)
print('    0.90:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'two-sided', confidence=0.95)
print('    0.95:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'two-sided', confidence=0.99)
print('    0.99:   pvalue =', res['p-val'][0], '  power =', res['power'][0])

print('E(X) > E(Y)')
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'greater', confidence=0.9)
print('    0.90:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'greater', confidence=0.95)
print('    0.95:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'greater', confidence=0.99)
print('    0.99:   pvalue =', res['p-val'][0], '  power =', res['power'][0])

print('E(X) < E(Y)')
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'less', confidence=0.9)
print('    0.90:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'less', confidence=0.95)
print('    0.95:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
res = ttest(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, paired = False, alternative = 'less', confidence=0.99)
print('    0.99:   pvalue =', res['p-val'][0], '  power =', res['power'][0])
mean(X) = 43.452991452991455 vs mean(Y) = 41.62068965517241
E(X) != E(Y)
    0.90:   pvalue = 0.2969828479314967   power = 0.1803430386782663
    0.95:   pvalue = 0.2969828479314967   power = 0.1803430386782663
    0.99:   pvalue = 0.2969828479314967   power = 0.1803430386782663
E(X) > E(Y)
    0.90:   pvalue = 0.14849142396574835   power = 0.2734624976802309
    0.95:   pvalue = 0.14849142396574835   power = 0.2734624976802309
    0.99:   pvalue = 0.14849142396574835   power = 0.2734624976802309
E(X) < E(Y)
    0.90:   pvalue = 0.8515085760342517   power = 0.0036011101618038895
    0.95:   pvalue = 0.8515085760342517   power = 0.0036011101618038895
    0.99:   pvalue = 0.8515085760342517   power = 0.0036011101618038895

R:

In [174]:
%%R

cat("mean(X) =", mean(df1[df1$gender %like% 'M', "age"]),
    "vs mean(Y) =", mean(mean(df1[df1$gender %like% 'F', "age"])))

cat("\nE(X) != E(Y)\n")
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "two.sided", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "two.sided", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "two.sided", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')

cat("\nE(X) > E(Y)\n")
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "greater", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "greater", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "greater", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')

cat("\nE(X) < E(Y)\n")
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "less", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "less", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = t.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], paired = FALSE, alternative = "less", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')
mean(X) = 43.45299 vs mean(Y) = 41.62069
E(X) != E(Y)
    0.90:  pvalue = 0.2969828 
    0.95:  pvalue = 0.2969828 
    0.99:  pvalue = 0.2969828 

E(X) > E(Y)
    0.90:  pvalue = 0.1484914 
    0.95:  pvalue = 0.1484914 
    0.99:  pvalue = 0.1484914 

E(X) < E(Y)
    0.90:  pvalue = 0.8515086 
    0.95:  pvalue = 0.8515086 
    0.99:  pvalue = 0.8515086 

Средний возраст авторов-мужчин выше среднего возраста авторов-женщин.

Критерий Уилкоксона - Манна - Уитни¶

Даны независимые выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m.$ $$H_0: \mathbb{E}X = \mathbb{E}Y,$$ т.е. матожидания выборок равны.

Альтернативные гипотезы: $$H_1: \mathbb{E}X \neq \mathbb{E}Y \quad (\text{двухсторонняя гипотеза})$$ $$H_1': \mathbb{E}X > \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$ $$H_1'': \mathbb{E}X < \mathbb{E}Y \quad (\text{односторонняя гипотеза})$$

Рассмотрим, как отличаются средние рейтинги комедий и драм (считаем выборки независимыми, так как у них всего 2 пересечения).

Python:

In [175]:
print('mean(X) =', statistics.mean(comedy_df.average_rating),
      'vs mean(Y) =', statistics.mean(drama_df.average_rating))

_, pvalue = mannwhitneyu(comedy_df.average_rating, drama_df.average_rating, alternative = 'two-sided')
print('E(X) != E(Y):   pvalue =', pvalue)

_, pvalue = mannwhitneyu(comedy_df.average_rating, drama_df.average_rating, alternative = 'greater')
print('E(X) > E(Y):   pvalue =', pvalue)

_, pvalue = mannwhitneyu(comedy_df.average_rating, drama_df.average_rating, alternative = 'less')
print('E(X) < E(Y):   pvalue =', pvalue)
mean(X) = 4.024275862068966 vs mean(Y) = 4.029975369458128
E(X) != E(Y):   pvalue = 0.38714692405450324
E(X) > E(Y):   pvalue = 0.8065933038563693
E(X) < E(Y):   pvalue = 0.19357346202725162

R:

In [176]:
%%R

cat("E(X) =", mean(comedy_df$average_rating),
    "vs E(Y) =", mean(mean(drama_df$average_rating)))

cat("\nE(X) == E(Y)\n")
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "two.sided", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "two.sided", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "two.sided", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')

cat("\nE(X) < E(Y)\n")
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "greater", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "greater", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "greater", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')

cat("\nE(X) > E(Y)\n")
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "less", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "less", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = wilcox.test(comedy_df$average_rating, drama_df$average_rating, paired = FALSE, alternative = "less", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')
E(X) = 4.024276 vs E(Y) = 4.029975
E(X) == E(Y)
    0.90:  pvalue = 0.3871469 
    0.95:  pvalue = 0.3871469 
    0.99:  pvalue = 0.3871469 

E(X) < E(Y)
    0.90:  pvalue = 0.8065933 
    0.95:  pvalue = 0.8065933 
    0.99:  pvalue = 0.8065933 

E(X) > E(Y)
    0.90:  pvalue = 0.1935735 
    0.95:  pvalue = 0.1935735 
    0.99:  pvalue = 0.1935735 

Средний рейтинг комедий меньше среднего рейтинга драм (более "серьёзные" произведения оцениваются выше).

Проверка гипотез об однородности дисперсий¶

Критерий Фишера¶

Даны выборки $X_1, \dots, X_n \sim N(\mu_1, \sigma_1)$ и $Y_1, \dots, Y_m \sim N(\mu_2, \sigma_2).$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.

Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$ $$H_1': DX > DY \quad (\text{односторонняя гипотеза})$$ $$H_1'': DX < DY \quad (\text{односторонняя гипотеза})$$

Рассмотрим, как отличаются разброс возрастов авторов-женщин и авторов-мужчин.

Python:


На Python нет реализации критерия Фишера, поэтому сделаем её сами.

In [177]:
def f_test(x, y, alternative='two_sided'):
    df1 = len(x) - 1
    df2 = len(y) - 1
    f = x.var() / y.var()
    if alternative == 'greater':
        p_val = 1.0 - stats.f.cdf(f, df1, df2)
    elif alternative == 'less':
        p_val = stats.f.cdf(f, df1, df2)
    else:
        p_val = 2.0 * (1.0 - stats.f.cdf(f, df1, df2))
    print('    pvalue = ', p_val)


print('var(X) =', statistics.variance(df1[df1.gender == 'M'].age),
      'vs var(Y) =', statistics.variance(df1[df1.gender == 'F'].age))

print('D(X) != D(Y):', end='')
f_test(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, alternative='two-sided')

print('D(X) > D(Y):', end='')
f_test(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, alternative='greater')

print('D(X) < D(Y):', end='')
f_test(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age, alternative='less')
var(X) = 119.11199528440908 vs var(Y) = 118.34482758620689
D(X) != D(Y):    pvalue =  0.9983269033717543
D(X) > D(Y):    pvalue =  0.49916345168587717
D(X) < D(Y):    pvalue =  0.5008365483141228

R:

In [178]:
%%R

cat("var(X) =", var(df1[df1$gender %like% 'M', "age"]),
    "vs var(Y) =", var(df1[df1$gender %like% 'F', "age"]))

cat("\nD(X) == D(Y)\n")
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "two.sided", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "two.sided", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "two.sided", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')

cat("\nD(X) > D(Y)\n")
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "greater", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "greater", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "greater", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')

cat("\nD(X) < D(Y)\n")
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "less", conf.level = 0.9)
cat("    0.90:  pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "less", conf.level = 0.95)
cat("    0.95:  pvalue =", res$p.value, '\n')
res = var.test(df1[df1$gender %like% 'M', "age"], df1[df1$gender %like% 'F', "age"], alternative = "less", conf.level = 0.99)
cat("    0.99:  pvalue =", res$p.value, '\n')
var(X) = 119.112 vs var(Y) = 118.3448
D(X) == D(Y)
    0.90:  pvalue = 0.9983269 
    0.95:  pvalue = 0.9983269 
    0.99:  pvalue = 0.9983269 

D(X) > D(Y)
    0.90:  pvalue = 0.4991635 
    0.95:  pvalue = 0.4991635 
    0.99:  pvalue = 0.4991635 

D(X) < D(Y)
    0.90:  pvalue = 0.5008365 
    0.95:  pvalue = 0.5008365 
    0.99:  pvalue = 0.5008365 

Дисперсия возрастов у мужчин и женщин равны.

Критерий Левене¶

Даны выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m.$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.

Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$

Рассмотрим, как отличаются дисперсии оценок комедий и драм.

Python:

In [179]:
print('var(X) =', np.var(comedy_df.average_rating),
      'vs var(Y) =', np.var(drama_df.average_rating))

_, pvalue = levene(comedy_df.average_rating, drama_df.average_rating)
print('D(X) != D(Y):   pvalue =', pvalue)
var(X) = 0.07031688941736028 vs var(Y) = 0.08652389101895216
D(X) != D(Y):   pvalue = 0.01462297954409574

R:

In [180]:
%%R

xor_copy <- df2[xor(df2$genres %like% "Comedy", df2$genres %like% "Drama"), ]
xor_copy[xor_copy$genres %like% "Comedy", "is_comedy"] <- 1
xor_copy[xor_copy$genres %like% "Drama", "is_comedy"] <- 0

and_copy <- df2[rep(which(df2$genres %like% "Comedy" & df2$genres %like% "Drama"), 2), ]
and_copy[c(1, 2), "is_comedy"] <- 1
and_copy[c(3, 4), "is_comedy"] <- 0

res_copy <- rbind(xor_copy, and_copy)
In [181]:
%%R

cat("var(X) =", var(comedy_df$average_rating),
      "vs var(Y) =", var(drama_df$average_rating), '\n\n')

leveneTest(res_copy$average_rating, res_copy$is_comedy)
var(X) = 0.0708052 vs var(Y) = 0.08673753 

Levene's Test for Homogeneity of Variance (center = median)
       Df F value  Pr(>F)  
group   1  5.9993 0.01462 *
      549                  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Дисперсии оценок драм и комедий не равны.

Критерий Бартлетта¶

Даны выборки $X_1, \dots, X_n \sim N(\mu_1, \sigma_1)$ и $Y_1, \dots, Y_m \sim N(\mu_2, \sigma_3).$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.

Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$ $$H_1': DX > DY \quad (\text{односторонняя гипотеза})$$ $$H_1'': DX < DY \quad (\text{односторонняя гипотеза})$$

Рассмотрим, как отличаются отклонения возрастов авторов-женщин и авторов-мужчин.

Python:

In [182]:
print('var(X) =', np.var(df1[df1.gender == 'M'].age),
      'vs var(Y) =', np.var(df1[df1.gender == 'F'].age))

_, pvalue = bartlett(df1[df1.gender == 'M'].age, df1[df1.gender == 'F'].age)

print('D(X) != D(Y):   pvalue =', pvalue)
var(X) = 118.09394404266199 vs var(Y) = 116.30439952437574
D(X) != D(Y):   pvalue = 0.9775499143207277

R:

In [183]:
%%R

cat("var(X) =", var(df1[df1$gender %like% "F", "age"]),
      "vs var(Y) =", var(df1[df1$gender %like% "M", "age"]))

res <- bartlett.test(df1$age ~ df1$gender)
cat("\nD(X) != D(Y):    pvalue =", res$p.value)
var(X) = 118.3448 vs var(Y) = 119.112
D(X) != D(Y):    pvalue = 0.9775499

Дисперсии возрастов у мужчин и женщин одинаковые.

Критерий Флигнера - Килина¶

Даны выборки $X_1, \dots, X_n$ и $Y_1, \dots, Y_m.$ $$H_0: DX = DY,$$ т.е. дисперсии выборок равны.

Альтернативные гипотезы: $$H_1: DX \neq DY \quad (\text{двухсторонняя гипотеза})$$ $$H_1': DX > DY \quad (\text{односторонняя гипотеза})$$ $$H_1'': DX < DY \quad (\text{односторонняя гипотеза})$$

Рассмотрим, как отличаются отклонения оценок комедий и драм.

Python:

In [184]:
print('var(X) =', np.var(comedy_df.average_rating),
      'vs var(Y) =', np.var(drama_df.average_rating))

_, pvalue = fligner(comedy_df.average_rating, drama_df.average_rating)
print('D(X) != D(Y):   pvalue =', pvalue)
var(X) = 0.07031688941736028 vs var(Y) = 0.08652389101895216
D(X) != D(Y):   pvalue = 0.012561899939530143

R:

In [185]:
%%R

cat("var(X) =", var(comedy_df$average_rating),
      "vs var(Y) =", var(drama_df$average_rating), '\n\n')

res <- fligner.test(res_copy$average_rating, res_copy$is_comedy)
cat("\nD(X) != D(Y):    pvalue =", res$p.value)
var(X) = 0.0708052 vs var(Y) = 0.08673753 


D(X) != D(Y):    pvalue = 0.0125619

Дисперсии оценок драм и комедий не равны.

8. Исследование корреляционных взаимосвязей¶

Рассмотрим, как коррелируют количественные данные в обоих датасетах.

Коэффициент корреляции Пирсона¶

Python:

In [186]:
print("df1:")
print("birth_year x age:")
print("    ", pearsonr(df1.birth_year, df1.age))
print("birth_year x first_published:")
print("    ", pearsonr(df1.birth_year, df1.first_published))
print("birth_year x sales:")
print("    ", pearsonr(df1.birth_year, df1.sales))
print("age x first_published:")
print("    ", pearsonr(df1.age, df1.first_published))
print("age x sales:")
print("    ", pearsonr(df1.age, df1.sales))
print("first_published x sales:")
print("    ", pearsonr(df1.first_published, df1.sales))

print("df2:")
print("average_rating x number_of_ratings:")
print("    ", pearsonr(df2.average_rating, df2.number_of_ratings))
df1:
birth_year x age:
     PearsonRResult(statistic=-0.16111193049698522, pvalue=0.033173980197517605)
birth_year x first_published:
     PearsonRResult(statistic=0.9858049406054875, pvalue=5.245782166441118e-136)
birth_year x sales:
     PearsonRResult(statistic=-0.1273678745473681, pvalue=0.093017306705826)
age x first_published:
     PearsonRResult(statistic=0.006876377241550452, pvalue=0.928036833440516)
age x sales:
     PearsonRResult(statistic=0.02408953122716002, pvalue=0.7516710557994336)
first_published x sales:
     PearsonRResult(statistic=-0.12495272423683916, pvalue=0.0994374716756319)
df2:
average_rating x number_of_ratings:
     PearsonRResult(statistic=0.028439105170956217, pvalue=0.004453284358105489)

R:

In [187]:
%%R

cat("df1:\n")
cat("birth_year x age:\n    ")
print(cor.test(df1$birth_year, df1$age, method="pearson"))
cat("birth_year x first_published:\n    ")
print(cor.test(df1$birth_year, df1$first_published, method="pearson"))
cat("birth_year x sales:\n    ")
print(cor.test(df1$birth_year, df1$sales, method="pearson"))
cat("age x first_published:\n    ")
print(cor.test(df1$age, df1$first_published, method="pearson"))
cat("age x sales:\n    ")
print(cor.test(df1$age, df1$sales, method="pearson"))
cat("first_published x sales:\n    ")
print(cor.test(df1$first_published, df1$sales, method="pearson"))

cat("\ndf2:\n")
cat("average_rating x number_of_ratings:\n    ")
print(cor.test(df2$average_rating, df2$number_of_ratings, method="pearson"))
df1:
birth_year x age:
    
	Pearson's product-moment correlation

data:  df1$birth_year and df1$age
t = -2.1471, df = 173, p-value = 0.03317
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.30223176 -0.01308145
sample estimates:
       cor 
-0.1611119 

birth_year x first_published:
    
	Pearson's product-moment correlation

data:  df1$birth_year and df1$first_published
t = 77.228, df = 173, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.9809071 0.9894530
sample estimates:
      cor 
0.9858049 

birth_year x sales:
    
	Pearson's product-moment correlation

data:  df1$birth_year and df1$sales
t = -1.689, df = 173, p-value = 0.09302
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.27059828  0.02137919
sample estimates:
       cor 
-0.1273679 

age x first_published:
    
	Pearson's product-moment correlation

data:  df1$age and df1$first_published
t = 0.090447, df = 173, p-value = 0.928
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1416112  0.1550613
sample estimates:
        cor 
0.006876377 

age x sales:
    
	Pearson's product-moment correlation

data:  df1$age and df1$sales
t = 0.31694, df = 173, p-value = 0.7517
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.1246992  0.1718187
sample estimates:
       cor 
0.02408953 

first_published x sales:
    
	Pearson's product-moment correlation

data:  df1$first_published and df1$sales
t = -1.6565, df = 173, p-value = 0.09944
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.26832226  0.02383215
sample estimates:
       cor 
-0.1249527 


df2:
average_rating x number_of_ratings:
    
	Pearson's product-moment correlation

data:  df2$average_rating and df2$number_of_ratings
t = 2.8448, df = 9998, p-value = 0.004453
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.008843965 0.048012413
sample estimates:
       cor 
0.02843911 

Коэффициент корреляции Спирмена¶

Python:

In [188]:
print("df1:")
print("birth_year x age:")
print("    ", spearmanr(df1.birth_year, df1.age))
print("birth_year x first_published:")
print("    ", spearmanr(df1.birth_year, df1.first_published))
print("birth_year x sales:")
print("    ", spearmanr(df1.birth_year, df1.sales))
print("age x first_published:")
print("    ", spearmanr(df1.age, df1.first_published))
print("age x sales:")
print("    ", spearmanr(df1.age, df1.sales))
print("first_published x sales:")
print("    ", spearmanr(df1.first_published, df1.sales))

print("df2:")
print("average_rating x number_of_ratings:")
print("    ", spearmanr(df2.average_rating, df2.number_of_ratings))
df1:
birth_year x age:
     SignificanceResult(statistic=-0.22161108578341485, pvalue=0.0032049243336212097)
birth_year x first_published:
     SignificanceResult(statistic=0.9286580663365899, pvalue=2.0088541923201725e-76)
birth_year x sales:
     SignificanceResult(statistic=-0.16708646133214528, pvalue=0.02710059167888117)
age x first_published:
     SignificanceResult(statistic=0.11232645210830353, pvalue=0.1388790585617484)
age x sales:
     SignificanceResult(statistic=0.025341821629909604, pvalue=0.7392159748274548)
first_published x sales:
     SignificanceResult(statistic=-0.1923388833606713, pvalue=0.010771516916827907)
df2:
average_rating x number_of_ratings:
     SignificanceResult(statistic=-0.15452971557803954, pvalue=1.733872942153623e-54)

R:

In [189]:
%%R

cat("df1:\n")
cat("birth_year x age:\n    ")
print(cor.test(df1$birth_year, df1$age, method="spearman"))
cat("birth_year x first_published:\n    ")
print(cor.test(df1$birth_year, df1$first_published, method="spearman"))
cat("birth_year x sales:\n    ")
print(cor.test(df1$birth_year, df1$sales, method="spearman"))
cat("age x first_published:\n    ")
print(cor.test(df1$age, df1$first_published, method="spearman"))
cat("age x sales:\n    ")
print(cor.test(df1$age, df1$sales, method="spearman"))
cat("first_published x sales:\n    ")
print(cor.test(df1$first_published, df1$sales, method="spearman"))

cat("\ndf2:\n")
cat("average_rating x number_of_ratings:\n    ")
print(cor.test(df2$average_rating, df2$number_of_ratings, method="spearman"))
df1:
birth_year x age:
    
	Spearman's rank correlation rho

data:  df1$birth_year and df1$age
S = 1091143, p-value = 0.003205
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.2216111 

birth_year x first_published:
    
	Spearman's rank correlation rho

data:  df1$birth_year and df1$first_published
S = 63723, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.9286581 

birth_year x sales:
    
	Spearman's rank correlation rho

data:  df1$birth_year and df1$sales
S = 1042442, p-value = 0.0271
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.1670865 

age x first_published:
    
	Spearman's rank correlation rho

data:  df1$age and df1$first_published
S = 792870, p-value = 0.1389
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.1123265 

age x sales:
    
	Spearman's rank correlation rho

data:  df1$age and df1$sales
S = 870565, p-value = 0.7392
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
0.02534182 

first_published x sales:
    
	Spearman's rank correlation rho

data:  df1$first_published and df1$sales
S = 1064997, p-value = 0.01077
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.1923389 


df2:
average_rating x number_of_ratings:
    
	Spearman's rank correlation rho

data:  df2$average_rating and df2$number_of_ratings
S = 1.9242e+11, p-value < 2.2e-16
alternative hypothesis: true rho is not equal to 0
sample estimates:
       rho 
-0.1545297 

Коэффициент корреляции Кендалла¶

Python:

In [190]:
print("df1:")
print("birth_year x age:")
print("    ", kendalltau(df1.birth_year, df1.age))
print("birth_year x first_published:")
print("    ", kendalltau(df1.birth_year, df1.first_published))
print("birth_year x sales:")
print("    ", kendalltau(df1.birth_year, df1.sales))
print("age x first_published:")
print("    ", kendalltau(df1.age, df1.first_published))
print("age x sales:")
print("    ", kendalltau(df1.age, df1.sales))
print("first_published x sales:")
print("    ", kendalltau(df1.first_published, df1.sales))

print("df2:")
print("average_rating x number_of_ratings:")
print("    ", kendalltau(df2.average_rating, df2.number_of_ratings))
df1:
birth_year x age:
     SignificanceResult(statistic=-0.1487202050324284, pvalue=0.004121070172677824)
birth_year x first_published:
     SignificanceResult(statistic=0.7829564425934482, pvalue=1.949443685647089e-52)
birth_year x sales:
     SignificanceResult(statistic=-0.11886744168201617, pvalue=0.02413018270093875)
age x first_published:
     SignificanceResult(statistic=0.07978447872205262, pvalue=0.12306653247813049)
age x sales:
     SignificanceResult(statistic=0.011638379456385595, pvalue=0.8264558116845667)
first_published x sales:
     SignificanceResult(statistic=-0.1350582730200036, pvalue=0.010250945964161366)
df2:
average_rating x number_of_ratings:
     SignificanceResult(statistic=-0.1030334892459667, pvalue=2.526856184366521e-53)

R:

In [191]:
%%R

cat("df1:\n")
cat("birth_year x age:\n    ")
print(cor.test(df1$birth_year, df1$age, method="kendall"))
cat("birth_year x first_published:\n    ")
print(cor.test(df1$birth_year, df1$first_published, method="kendall"))
cat("birth_year x sales:\n    ")
print(cor.test(df1$birth_year, df1$sales, method="kendall"))
cat("age x first_published:\n    ")
print(cor.test(df1$age, df1$first_published, method="kendall"))
cat("age x sales:\n    ")
print(cor.test(df1$age, df1$sales, method="kendall"))
cat("first_published x sales:\n    ")
print(cor.test(df1$first_published, df1$sales, method="kendall"))

cat("\ndf2:\n")
cat("average_rating x number_of_ratings:\n    ")
print(cor.test(df2$average_rating, df2$number_of_ratings, method="kendall"))
df1:
birth_year x age:
    
	Kendall's rank correlation tau

data:  df1$birth_year and df1$age
z = -2.8687, p-value = 0.004121
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
-0.1487202 

birth_year x first_published:
    
	Kendall's rank correlation tau

data:  df1$birth_year and df1$first_published
z = 15.239, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
      tau 
0.7829564 

birth_year x sales:
    
	Kendall's rank correlation tau

data:  df1$birth_year and df1$sales
z = -2.2551, p-value = 0.02413
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
-0.1188674 

age x first_published:
    
	Kendall's rank correlation tau

data:  df1$age and df1$first_published
z = 1.542, p-value = 0.1231
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
0.07978448 

age x sales:
    
	Kendall's rank correlation tau

data:  df1$age and df1$sales
z = 0.21925, p-value = 0.8265
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
0.01163838 

first_published x sales:
    
	Kendall's rank correlation tau

data:  df1$first_published and df1$sales
z = -2.5672, p-value = 0.01025
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
-0.1350583 


df2:
average_rating x number_of_ratings:
    
	Kendall's rank correlation tau

data:  df2$average_rating and df2$number_of_ratings
z = -15.372, p-value < 2.2e-16
alternative hypothesis: true tau is not equal to 0
sample estimates:
       tau 
-0.1030335 

Видна высокая корреляция между birth_year и first_published в первом датасете.

9. Использование методов $\chi^2$, точного теста Фишера, теста МакНемара, Кохрана - Мантеля - Хензеля¶

Тест $\chi^2$¶

Тест $\chi^2$ используется, для проверки наличия связи между двумя категориальными признаками ($X$ и $Y$), первый из которых принимает $k$ значений, а второй - $l$. Соответствующие частоты заносят в таблицу сопряжённости:

$x_1$ $x_2$ ... $x_l$
$y_1$ $f_{11}$ $f_{12}$ ... $f_{1l}$
$y_2$ $f_{21}$ $f_{22}$ ... $f_{2l}$
... ...
$y_k$ $f_{k1}$ $f_{k2}$ ... $f_{kl}$
$$H_0: \mathbb{P}\{\ X=x_i,\ Y=y_j\ \} = \mathbb{P}\{\ X=x_i\ \} * \mathbb{P}\{\ Y=y_j\ \}, \quad 1\leq i\leq k, 1\leq j\leq l,$$

т.е. переменные $X$ и $Y$ независимы.

Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.

Проверим, зависит ли рейтинг от жанра. Для этого разделим шкалу оценивания (0 - 5) на 5 равных частей и cоставим таблицу сопряжённости.

Python:

In [230]:
# делим шкалу оценивания
df2['rating_cut'] = pd.cut(df2['average_rating'], bins=[-0.1, 1, 2, 3, 4, 5], labels=[1, 2, 3, 4, 5])
# отделяем жанры друг от друга
df2_explode = df2.explode('genres')
# составляем таблицу сопряжённости
cross_df = pd.crosstab(df2_explode.genres, df2_explode.rating_cut)
print('Количество строк:', cross_df.shape[0])
cross_df.head(5)
Количество строк: 617
Out[230]:
rating_cut 2 3 4 5
genres
12th Century 0 0 0 1
15th Century 0 0 1 1
16th Century 0 0 6 2
17th Century 0 0 4 3
18th Century 0 0 12 4
In [193]:
chi2_contingency(cross_df)
Out[193]:
Chi2ContingencyResult(statistic=5118.839776700974, pvalue=2.868874482754211e-304, dof=1848, expected_freq=array([[6.90095406e-05, 9.14376413e-04, 4.34242534e-01, 5.64774080e-01],
       [1.38019081e-04, 1.82875283e-03, 8.68485068e-01, 1.12954816e+00],
       [5.52076325e-04, 7.31501130e-03, 3.47394027e+00, 4.51819264e+00],
       ...,
       [2.76038162e-04, 3.65750565e-03, 1.73697014e+00, 2.25909632e+00],
       [6.90095406e-05, 9.14376413e-04, 4.34242534e-01, 5.64774080e-01],
       [2.89840070e-03, 3.84038093e-02, 1.82381864e+01, 2.37205114e+01]]))

R:

In [231]:
%%R

# делим шкалу оценивания
df2$rating_cut <- cut(df2$average_rating, breaks=c(-0.1, 1, 2, 3, 4, 5), labels=c(1, 2, 3, 4, 5))
# отделяем жанры друг от друга
df2_explode <- df2 %>% tidyr::unnest(genres)
# составляем таблицу сопряжённости
cross_df <- table(df2_explode$genres, df2_explode$rating_cut)
cat("Количество строк:", nrow(cross_df), "\n")
head(cross_df, n=5)
Количество строк: 617 
             
               1  2  3  4  5
  12thCentury  0  0  0  0  1
  15thCentury  0  0  0  1  1
  16thCentury  0  0  0  6  2
  17thCentury  0  0  0  4  3
  18thCentury  0  0  0 12  4
In [195]:
%%R

chisq.test(cross_df)
	Pearson's Chi-squared test

data:  cross_df
X-squared = NaN, df = 2464, p-value = NA

Рейтинг и жанр зависимы.

Точный тест Фишера¶

Точный тест Фишера используется, для проверки наличия связи между двумя категориальными признаками ($X$ и $Y$), каждый из которых принимает по $2$ значения:

$x_1$ $x_2$
$y_1$ $f_{11}$ $f_{12}$
$y_2$ $f_{21}$ $f_{22}$
$$H_0: \mathbb{P}\{\ X=x_i,\ Y=y_j\ \} = \mathbb{P}\{\ X=x_i\ \} * \mathbb{P}\{\ Y=y_j\ \}, \quad 1\leq i\leq 2,\quad 1\leq j\leq 2,$$

т.е. переменные $X$ и $Y$ независимы.

Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.

Фишера можно применять на меньших выборках данных, чем критерий $\chi^2$. Воспользуемся этим, чтобы найти пересечение используемых датасетов и проверить, существует ли связь между оценкой, которую получило произведение (больше или меньше медианы), и его продажами (больше или меньше медианы).

Python:

In [196]:
ultra_df = pd.merge(df1, df2, how='inner', on=['title']).drop(['author_y', 'genres_x', 'rating_cut'], axis= 1).rename(columns = {'author_x':'author', 'genres_y':'genres'})
# делим шкалы оценивания
ultra_df['rating_cut'] = pd.cut(ultra_df['average_rating'],
                                bins=[min(ultra_df.average_rating) - 1, statistics.median(ultra_df.average_rating),
                                      max(ultra_df.average_rating)],
                                labels=['rating: less', 'rating: greater'])
ultra_df['sales_cut'] = pd.cut(ultra_df['sales'],
                                bins=[min(ultra_df.sales) - 1, statistics.median(ultra_df.sales), max(ultra_df.sales)],
                               labels=['sales: less', 'sales: greater'])
# составляем таблицу сопряжённости
cross_df = pd.crosstab(ultra_df.sales_cut, ultra_df.rating_cut)
cross_df
Out[196]:
rating_cut rating: less rating: greater
sales_cut
sales: less 18 19
sales: greater 18 14
In [197]:
fisher_exact(cross_df)
Out[197]:
SignificanceResult(statistic=0.7368421052631579, pvalue=0.6308158496561179)

R:

In [198]:
%%R

ultra_df <- subset(merge(df1, df2, by = "title"), select = -c(genres.x, author.y, rating_cut))
names(ultra_df)[names(ultra_df) == "author.x"] <- "author"
names(ultra_df)[names(ultra_df) == "genres.y"] <- "genres"
# делим шкалы оценивания
ultra_df$rating_cut <- cut(ultra_df$average_rating,
                           breaks=c(min(ultra_df$average_rating) - 1, median(ultra_df$average_rating),
                                    max(ultra_df$average_rating)),
                           labels=c("rating: less", "rating: greater"))
ultra_df$sales_cut <- cut(ultra_df$sales,
                          breaks=c(min(ultra_df$sales) - 1, median(ultra_df$sales),
                                   max(ultra_df$sales)),
                          labels=c("sales: less", "sales: greater"))
# составляем таблицу сопряжённости
cross_df <- table(ultra_df$sales_cut, ultra_df$rating_cut)
cross_df
                
                 rating: less rating: greater
  sales: less              18              19
  sales: greater           18              14
In [199]:
%%R

fisher.test(cross_df)
	Fisher's Exact Test for Count Data

data:  cross_df
p-value = 0.6308
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.2558414 2.1118582
sample estimates:
odds ratio 
 0.7401323 

Переменые независимы, связи между рейтингом и количеством проданных копий нет.

Тест МакНемара¶

Когда независимости наблюдений нет (признак проверяется на одних и тех же субъектах), для анализа таблиц напряжённости 2x2 используется критерий МакНемара.

Рассмотрим $n$ субъектов, для каждого из которых проведено 2 теста с 2 исходами ("+" и "-"). Результаты тестов занесены в таблицу сопряжённости:

$\text{Тест 1}$ $\text{ВСЕГО}$
$+$
$-$
$\text{Тест 2}$ $+$ $a$ $b$ $a + b$
$-$ $c$ $d$ $c + d$
$\text{ВСЕГО}$ $a + c$ $b + d$ $a + b + c + d = n$
$$H_0: \quad p_a + p_b = p_a + p_c, \quad p_c + p_d = p_b + p_d, \quad(\Rightarrow \ \ p_b = p_c)$$

т.е. маргинальные распределения для всех исходов совпадают.

Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.

Рассмотрим, изменилось ли со временем влияние пола автора на жанр произведения. Для этого найдём количество жанров, в которых большинство произведений написаны мужчиной/женщиной до/после 1990.

Python:

In [200]:
df1_explode = df1.explode('genres')

# до 1990
fem_genres_freq = df1_explode[(df1_explode.gender == 'F') & (df1_explode.first_published < 1990)].genres.value_counts().reset_index()
fem_genres_freq.columns = ['genre', 'fem_freq_before']

male_genres_freq = df1_explode[(df1_explode.gender == 'M') & (df1_explode.first_published < 1990)].genres.value_counts().reset_index()
male_genres_freq.columns = ['genre', 'male_freq_before']

genres_freq_before = pd.merge(fem_genres_freq, male_genres_freq, 'outer')

# после 1990
fem_genres_freq = df1_explode[(df1_explode.gender == 'F') & (df1_explode.first_published >= 1990)].genres.value_counts().reset_index()
fem_genres_freq.columns = ['genre', 'fem_freq_after']

male_genres_freq = df1_explode[(df1_explode.gender == 'M') & (df1_explode.first_published >= 1990)].genres.value_counts().reset_index()
male_genres_freq.columns = ['genre', 'male_freq_after']

genres_freq_after = pd.merge(fem_genres_freq, male_genres_freq, 'outer')

# объединяем и суммируем
genres_freq = pd.merge(genres_freq_before, genres_freq_after, 'outer').fillna(0)

genres_freq['fem_win_before'] = np.where(genres_freq.fem_freq_before >= genres_freq.male_freq_before, 1, 0)
genres_freq['fem_win_after'] = np.where(genres_freq.fem_freq_after >= genres_freq.male_freq_after, 1, 0)

genres_freq['win_win'] = np.where((genres_freq.fem_win_before == 1) & (genres_freq.fem_win_after == 1), 1, 0)
genres_freq['win_lose'] = np.where((genres_freq.fem_win_before == 1) & (genres_freq.fem_win_after == 0), 1, 0)
genres_freq['lose_win'] = np.where((genres_freq.fem_win_before == 0) & (genres_freq.fem_win_after == 1), 1, 0)
genres_freq['lose_lose'] = np.where((genres_freq.fem_win_before == 0) & (genres_freq.fem_win_after == 0), 1, 0)

genres_freq.head(n=5)
Out[200]:
genre fem_freq_before male_freq_before fem_freq_after male_freq_after fem_win_before fem_win_after win_win win_lose lose_win lose_lose
0 Fiction 23.0 67.0 23.0 20.0 0 1 0 0 1 0
1 Classics 22.0 68.0 4.0 10.0 0 0 0 0 0 1
2 Childrens 10.0 17.0 8.0 1.0 0 1 0 0 1 0
3 Historical Fiction 10.0 25.0 2.0 9.0 0 0 0 0 0 1
4 Romance 8.0 5.0 9.0 3.0 1 1 1 0 0 0
In [201]:
data = []
data.append([genres_freq['win_win'].sum(), genres_freq.shape[0] - genres_freq['win_lose'].sum()])
data.append([genres_freq['lose_win'].sum(), genres_freq.shape[0] - genres_freq['lose_lose'].sum()])
data
Out[201]:
[[41, 153], [79, 138]]
In [202]:
print(mcnemar(data, exact=False))
pvalue      1.6456406755618912e-06
statistic   22.969827586206897

R:

In [203]:
%%R

df1_explode <- df1 %>% tidyr::unnest(genres)

# до 1990
fem_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'F') & (df1_explode$first_published < 1990),]$genres))
names(fem_genres_freq) <- c("genre", "fem_freq_before")

male_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'M') & (df1_explode$first_published < 1990),]$genres))
names(male_genres_freq) <- c("genre", "male_freq_before")

genres_freq_before <- merge(fem_genres_freq, male_genres_freq, all = TRUE)

# после 1990
fem_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'F') & (df1_explode$first_published >= 1990),]$genres))
names(fem_genres_freq) <- c("genre", "fem_freq_after")

male_genres_freq <- as.data.frame(table(df1_explode[(df1_explode$gender %like% 'M') & (df1_explode$first_published >= 1990),]$genres))
names(male_genres_freq) <- c("genre", "male_freq_after")

genres_freq_after <- merge(fem_genres_freq, male_genres_freq, all = TRUE)

# объединяем и суммируем
genres_freq <- merge(genres_freq_before, genres_freq_after, all = TRUE)
genres_freq[is.na(genres_freq)] <- 0

genres_freq$fem_win_before <- ifelse(genres_freq$fem_freq_before >= genres_freq$male_freq_before, 1, 0)
genres_freq$fem_win_after <- ifelse(genres_freq$fem_freq_after >= genres_freq$male_freq_after, 1, 0)

genres_freq$win_win <- ifelse(genres_freq$fem_win_before == 1 & genres_freq$fem_win_after == 1, 1, 0)
genres_freq$win_lose <- ifelse(genres_freq$fem_win_before == 1 & genres_freq$fem_win_after == 0, 1, 0)
genres_freq$lose_win <- ifelse(genres_freq$fem_win_before == 0 & genres_freq$fem_win_after == 1, 1, 0)
genres_freq$lose_lose <- ifelse(genres_freq$fem_win_before == 0 & genres_freq$fem_win_after == 0, 1, 0)

head(genres_freq, n=5)
            genre fem_freq_before male_freq_before fem_freq_after
1           Adult               4                0              6
2       Adventure               1               15              6
3         Animals               4                4              1
4            Asia               1                2              1
5 AsianLiterature               1                1              0
  male_freq_after fem_win_before fem_win_after win_win win_lose lose_win
1               4              1             1       1        0        0
2               5              0             1       0        0        1
3               4              1             0       0        1        0
4               1              0             1       0        0        1
5               0              1             1       1        0        0
  lose_lose
1         0
2         0
3         0
4         0
5         0
In [204]:
%%R

data <- list()
data[[1]] <- c(sum(genres_freq$win_win), nrow(genres_freq) - sum(genres_freq$win_lose))
data[[2]] <- c(sum(genres_freq$lose_win), nrow(genres_freq) - sum(genres_freq$lose_lose))
In [205]:
%%R

mcnemar.test(matrix(unlist(data), ncol=2, byrow=TRUE))
	McNemar's Chi-squared test with continuity correction

data:  matrix(unlist(data), ncol = 2, byrow = TRUE)
McNemar's chi-squared = 22.97, df = 1, p-value = 1.646e-06

Со временем количество женров, в которых пишут преимущественно женщины, изменилось.

Тест Кохрана - Мантеля - Хензеля¶

Если наблюдения сгруппированы по $k$ стратам, то используется тест Кохрана-Мантеля-Хензеля. В таком случае рассматривают $k$ таблиц сопряжённости 2x2:

$x_1$ $x_2$ $\text{ВСЕГО}$
$y_1$ $a_i$ $b_i$ $a_i+b_i$
$y_2$ $c_i$ $d_i$ $c_i+d_i$
$\text{ВСЕГО}$ $a_i+c_i$ $b_i+d_i$ $t_i=a_i+b_i+c_i+d_i$
$$H_0: R = \frac{\sum^k_{i=1} p_{a_i} / p_{b_i}}{\sum^k_{i=1} p_{c_i} / p_{d_i}} = 1,$$

т.е. шансы равны (между признаками нет связи).

Если полученный с помощью критерия $pvalue$ меньше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.

Проверим, как связаны пол автора и жанр литературы (художестенная/нехудожественная). При этом поделим наблюдения на 3 части по возрасту автора: молодёжь (до 30 лет), взрослые (от 31 до 55 лет) и пожилые (старше 55 лет).

Python:

In [206]:
df1_explode['age_cut'] = pd.cut(df1_explode['age'],
                                bins=[min(df1_explode.age) - 1, 30,
                                      55, max(df1_explode.age)],
                                labels=['young', 'adult', 'senior'])
In [207]:
result = CMH(df1_explode[(df1_explode.genres == 'Nonfiction') |
            (df1_explode.genres == 'Fiction')],
             'gender', 'genres', stratifier='age_cut')
print(result)
        Cochran-Mantel-Haenszel Chi2 test

"gender" x "genres", stratified by "age_cut"

Cochran-Mantel-Haenszel M^2 = 0.2914, dof = 1, p-value = 0.5893

R:

In [208]:
%%R

df1_explode$age_cut <- cut(df1_explode$age,
                           breaks=c(min(df1_explode$age) - 1, 30,
                                    55, max(df1_explode$age)),
                           labels=c("young", "adult", "senior"))

mantelhaen_array <- array(
    c(nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "young", ]),
      nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "young", ]),
      nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "young", ]),
      nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "young", ]),

      nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "adult", ]),
      nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "adult", ]),
      nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "adult", ]),
      nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "adult", ]),

      nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "senior", ]),
      nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Fiction" & df1_explode$age_cut == "senior", ]),
      nrow(df1_explode[df1_explode$gender == "F" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "senior", ]),
      nrow(df1_explode[df1_explode$gender == "M" & df1_explode$genres == "Nonfiction" & df1_explode$age_cut == "senior", ])),
    dim = c(2, 2, 3),
    dimnames = list(
      gender = c("F", "M"),
      genres = c("Fiction", "Nonfiction"),
      age_cut = c("young", "adult", "senior")))
In [209]:
%%R

mantelhaen.test(mantelhaen_array)
	Mantel-Haenszel chi-squared test with continuity correction

data:  mantelhaen_array
Mantel-Haenszel X-squared = 0.11742, df = 1, p-value = 0.7318
alternative hypothesis: true common odds ratio is not equal to 1
95 percent confidence interval:
 0.5676904 2.6731700
sample estimates:
common odds ratio 
         1.231882 

Между полом автора и выбором жанра нет связи.

10. Проверка наличия мультиколлинеарности в данных¶

Корреляционная матрица¶

Простейший способ графически изобразить попарную корреляцию между величинами - корреляционная матрица. Построим такие матрицы для обоих датасетов.

Python:

In [210]:
data1 = df1[['birth_year', 'age', 'first_published', 'sales']].dropna()

fig = plt.figure(figsize=(7,6))

ax = sns.heatmap(data1.corr(),
            vmin=-1, vmax=1, cmap='RdBu',
            annot = True,
            annot_kws={"fontsize":9})
plt.xticks(fontsize=10)
plt.yticks(fontsize=10)
cbar = ax.collections[0].colorbar
cbar.ax.tick_params(labelsize=10)
plt.title('Корреляционная матрица: df1', fontsize=10, fontweight='bold')

plt.show()
In [211]:
data2 = df2[['average_rating', 'number_of_ratings']].dropna()

fig = plt.figure(figsize=(5,4))

ax = sns.heatmap(data2.corr(),
            vmin=-1, vmax=1, cmap='RdBu',
            annot = True,
            annot_kws={"fontsize":9})
plt.xticks(fontsize=10)
plt.yticks(fontsize=10)
cbar = ax.collections[0].colorbar
cbar.ax.tick_params(labelsize=10)
plt.title('Корреляционная матрица: df2', fontsize=10, fontweight='bold')

plt.show()

R:

In [212]:
%%R -w 600 -h 600

data1 <- dropNA(df1[c("birth_year", "age", "first_published", "sales")])

corrplot(cor(data1),
         method = "color",
         tl.col = "black", addCoef.col = "gray45")
In [213]:
%%R -w 500 -h 500

data2 <- dropNA(df2[c('average_rating', 'number_of_ratings')])

corrplot(cor(data2),
         method = "color",
         tl.col = "black", addCoef.col = "gray45")

В первом датасете, ожидаемо, нашлась высокая корреляция между годом первой публикации произведения и годом рождения его автора.

Фактор инфляции дисперсии¶

Одним из способов обнаружения мультиколлинеарности является использование коэффициента инфляции дисперсии (VIF). Найдём его для обоих датасетов.

Python:

In [214]:
print("df1:")

X = add_constant(data1)
ds = pd.Series([variance_inflation_factor(X.values, i) for i in range(X.shape[1])],
               index=X.columns)
print(ds)
df1:
const              981.311602
birth_year                inf
age                       inf
first_published           inf
sales                1.016504
dtype: float64
In [215]:
print("df2:")

X = add_constant(data2)
ds = pd.Series([variance_inflation_factor(X.values, i) for i in range(X.shape[1])],
               index=X.columns)
print(ds)
df2:
const                148.205952
average_rating         1.000809
number_of_ratings      1.000809
dtype: float64

R:

In [216]:
%%R

cat("df1:\n")

data1$const <- 1

# так как age и birth_year - коэффициенты-псевдонимы, используем только один из них за раз

mod <- lm(const ~ age + first_published + sales, data = data1)
cat("\n без birth_year:\n")
print(vif(mod))

mod <- lm(const ~ birth_year + first_published + sales, data = data1)
cat("\n без age:\n")
print(vif(mod))
df1:

 без birth_year:
            age first_published           sales 
       1.000680        1.015962        1.016504 

 без age:
     birth_year first_published           sales 
      35.497757       35.475782        1.016504 
In [217]:
%%R

cat("df2:\n")

data2$const <- 1

mod <- lm(const ~ average_rating + number_of_ratings, data = data2)
print(vif(mod))
df2:
   average_rating number_of_ratings 
         1.000809          1.000809 

VIF в первом датасете указывает на связь между birth_year, age и first_published (действительно, age = first_published - birth_year). Полученные с помощью VIF данные подтверждают то, что заметили с помощью матриц корреляции и коэффициентов корреляции.

11. Дисперсионный анализ¶

Чтобы исследовать влияние одной или нескольких переменных категориальных переменных (факторов) на одну зависимую количественную переменную (отклик), используется дисперсионный анализ (ANOVA).

Важно, чтобы дисперсии в группах, выделенных в связи с уровнями фактора, были равны и распределение было нормальным (или достаточно близким к нему).

Однофакторный дисперсионный анализ¶

Рассмотрим влияние фактора с $k$ уровнями. В соответствии с уровнями фактора выделяют $k$ подвыборок размерами $n_1, \dots, n_k$ соответственно: $ X_1 = (X_{11}, \dots, X_{1n_1});\ \dots;\ X_k = (X_{k1}, \dots, X_{kn_k}).$

$$H_0: \mathbb{E}(X_1) = \dots = \mathbb{E}(X_k),$$

т.е. средние выборок равны между собой (фактор не оказывает значимого влияния).

Если полученный с помощью критерия $pvalue$ больше 0.05, то нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.

Рассмотрим, зависит ли количество оценок от рейтинга.

Python:

In [218]:
group1 = df2.loc[df2.rating_cut == 1].number_of_ratings
group2 = df2.loc[df2.rating_cut == 2].number_of_ratings
group3 = df2.loc[df2.rating_cut == 3].number_of_ratings
group4 = df2.loc[df2.rating_cut == 4].number_of_ratings
group5 = df2.loc[df2.rating_cut == 5].number_of_ratings

f_oneway(group1, group2, group3, group4, group5)
Out[218]:
F_onewayResult(statistic=3.2063572589018983, pvalue=0.012197597739233967)

R:

In [219]:
%%R

summary(aov(number_of_ratings ~ rating_cut, df2))
              Df    Sum Sq  Mean Sq F value Pr(>F)  
rating_cut     4 1.500e+12 3.75e+11   3.206 0.0122 *
Residuals   9995 1.169e+15 1.17e+11                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Рейтинг книги влияет на количество оценок. Действительно, только произведения с маленьким количеством оценок получают низкий рейтинг:

In [220]:
print(max(df2[df2.rating_cut == 1].number_of_ratings))
print(max(df2[df2.rating_cut == 2].number_of_ratings))
print(max(df2[df2.rating_cut == 3].number_of_ratings))
print(max(df2[df2.rating_cut == 4].number_of_ratings))
print(max(df2[df2.rating_cut == 5].number_of_ratings))
0
2302
58521
6158643
9278135
In [221]:
%%R

cat(max(df2[df2$rating_cut == 1, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 2, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 3, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 4, "number_of_ratings"]), "\n")
cat(max(df2[df2$rating_cut == 5, "number_of_ratings"]), "\n")
0 
2302 
58521 
6158643 
9278135 

Многофакторный дисперсионный анализ¶

Опишем случай двухфакторного анализа (большее количество факторов аналогично).

Рассмотрим влияние факторов $A$ и $B$ с $k$ и $l$ уровнями. В соответствии с уровнями факторов выделяют $k*l$ подвыборок: $ X_{1,1}, \dots, X_{k,l}.$ В таком случае формулируют 3 нулевых гипотезы: $$H_0: \mathbb{E}(X_{1i}) = \dots = \mathbb{E}(X_{ki}), \forall i$$ т.е. средние под влиянием фактора $A$ равны (фактор $A$ не оказывает значимого влияния). $$H_0: \mathbb{E}(X_{i1}) = \dots = \mathbb{E}(X_{il}), \forall i$$ т.е. средние под влиянием фактора $B$ равны (фактор $B$ не оказывает значимого влияния).

$H_0:$ факторы $A$ и $B$ не взаимодействуют.

Если полученный с помощью критерия $pvalue$ больше 0.05, то соответствующая нулевая гипотеза отвергается, иначе - нет оснований отбросить гипотезу.

Проверим, влияют ли язык и пол и возраст автора на продажи.

Python:

In [222]:
df1['age_cut'] = pd.cut(df1['age'],
                        bins=[min(df1.age) - 1, 30,
                              55, max(df1.age)],
                        labels=['young', 'adult', 'senior'])
In [223]:
model = ols('sales ~ language + gender + age_cut + language:gender + language:age_cut + gender:age_cut + language:gender:age_cut', data=df1).fit()
sm.stats.anova_lm(model, typ=2)
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 15, but rank is 4
  warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 1, but rank is 0
  warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 2, but rank is 0
  warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 15, but rank is 4
  warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 30, but rank is 6
  warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 2, but rank is 0
  warnings.warn('covariance of constraints does not have full '
/usr/local/lib/python3.10/dist-packages/statsmodels/base/model.py:1888: ValueWarning: covariance of constraints does not have full rank. The number of constraints is 30, but rank is 16
  warnings.warn('covariance of constraints does not have full '
Out[223]:
sum_sq df F PR(>F)
language 318.481500 15.0 0.025895 0.998688
gender NaN 1.0 NaN NaN
age_cut NaN 2.0 NaN NaN
language:gender 7208.630522 15.0 0.586121 0.673192
language:age_cut 9493.975688 30.0 0.385969 0.887081
gender:age_cut NaN 2.0 NaN NaN
language:gender:age_cut 17166.811598 30.0 0.697901 0.792587
Residual 118069.326217 144.0 NaN NaN

R:

In [224]:
%%R

df1$age_cut <- cut(df1$age,
                   breaks=c(min(df1$age) - 1, 30,
                            55, max(df1$age)),
                   labels=c("young", "adult", "senior"))
In [225]:
%%R

summary(aov(sales ~ language * gender * age_cut, df1))
                  Df Sum Sq Mean Sq F value Pr(>F)
language          15   7179   478.6   0.584  0.884
gender             1     27    27.5   0.033  0.855
age_cut            2   2181  1090.4   1.330  0.268
language:gender    5   2733   546.6   0.667  0.649
language:age_cut   5   4838   967.6   1.180  0.322
gender:age_cut     2    299   149.6   0.183  0.833
Residuals        144 118069   819.9               

Вышеперечисленные факторы и их сочетания не влияют на продажи.

12. Регрессионные модели¶

Попробуем понять, как количество оценок влияет на рейтинг. Обозначим эти переменные как $y$ и $x$ соответственно. Сравним результаты регрессии методом наименьших квадратов следующих линейных по параметру моделей: линейной ($y=wx+b$), гиберболы ($y=w\frac{1}{x+1}+b$), корня ($y=w\sqrt{x}+b$) и логарифмической ($y=w\ln (x+1)+b$).

Python:

In [226]:
x = df2.number_of_ratings
y = df2.average_rating
x_no0 = df2[(df2.number_of_ratings != 0) & (df2.average_rating != 0)].number_of_ratings
y_no0 = df2[(df2.number_of_ratings != 0) & (df2.average_rating != 0)].average_rating

lin_model = sm.OLS(y, sm.add_constant(x)).fit()
print(lin_model.summary())

#X = np.column_stack((np.ones(x.shape[0]), x, x**2, x**3))
hyper_model = sm.OLS(y, sm.add_constant(1 / (x + 1))).fit()
print(hyper_model.summary())

sqrt_model = sm.OLS(y, sm.add_constant(x ** 0.5)).fit()
print(sqrt_model.summary())

log_model = sm.OLS(y, sm.add_constant(np.log(x + 1))).fit()
print(log_model.summary())
                            OLS Regression Results                            
==============================================================================
Dep. Variable:         average_rating   R-squared:                       0.001
Model:                            OLS   Adj. R-squared:                  0.001
Method:                 Least Squares   F-statistic:                     8.093
Date:                Mon, 11 Dec 2023   Prob (F-statistic):            0.00445
Time:                        12:34:39   Log-Likelihood:                -3259.3
No. Observations:               10000   AIC:                             6523.
Df Residuals:                    9998   BIC:                             6537.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
=====================================================================================
                        coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------------
const                 4.0660      0.003   1170.342      0.000       4.059       4.073
number_of_ratings  2.787e-08    9.8e-09      2.845      0.004    8.67e-09    4.71e-08
==============================================================================
Omnibus:                     4767.098   Durbin-Watson:                   1.932
Prob(Omnibus):                  0.000   Jarque-Bera (JB):           171661.269
Skew:                          -1.633   Prob(JB):                         0.00
Kurtosis:                      23.033   Cond. No.                     3.67e+05
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 3.67e+05. This might indicate that there are
strong multicollinearity or other numerical problems.
                            OLS Regression Results                            
==============================================================================
Dep. Variable:         average_rating   R-squared:                       0.002
Model:                            OLS   Adj. R-squared:                  0.002
Method:                 Least Squares   F-statistic:                     19.31
Date:                Mon, 11 Dec 2023   Prob (F-statistic):           1.12e-05
Time:                        12:34:39   Log-Likelihood:                -3253.7
No. Observations:               10000   AIC:                             6511.
Df Residuals:                    9998   BIC:                             6526.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
=====================================================================================
                        coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------------
const                 4.0648      0.003   1175.209      0.000       4.058       4.072
number_of_ratings     0.2544      0.058      4.395      0.000       0.141       0.368
==============================================================================
Omnibus:                     5701.374   Durbin-Watson:                   1.932
Prob(Omnibus):                  0.000   Jarque-Bera (JB):           272907.764
Skew:                          -2.042   Prob(JB):                         0.00
Kurtosis:                      28.264   Cond. No.                         17.3
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
                            OLS Regression Results                            
==============================================================================
Dep. Variable:         average_rating   R-squared:                       0.001
Model:                            OLS   Adj. R-squared:                  0.001
Method:                 Least Squares   F-statistic:                     6.808
Date:                Mon, 11 Dec 2023   Prob (F-statistic):            0.00909
Time:                        12:34:39   Log-Likelihood:                -3259.9
No. Observations:               10000   AIC:                             6524.
Df Residuals:                    9998   BIC:                             6538.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
=====================================================================================
                        coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------------
const                 4.0755      0.004    955.063      0.000       4.067       4.084
number_of_ratings  -3.65e-05    1.4e-05     -2.609      0.009   -6.39e-05   -9.08e-06
==============================================================================
Omnibus:                     4846.069   Durbin-Watson:                   1.930
Prob(Omnibus):                  0.000   Jarque-Bera (JB):           174716.871
Skew:                          -1.674   Prob(JB):                         0.00
Kurtosis:                      23.202   Cond. No.                         388.
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
                            OLS Regression Results                            
==============================================================================
Dep. Variable:         average_rating   R-squared:                       0.030
Model:                            OLS   Adj. R-squared:                  0.030
Method:                 Least Squares   F-statistic:                     310.9
Date:                Mon, 11 Dec 2023   Prob (F-statistic):           1.54e-68
Time:                        12:34:39   Log-Likelihood:                -3110.2
No. Observations:               10000   AIC:                             6224.
Df Residuals:                    9998   BIC:                             6239.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
=====================================================================================
                        coef    std err          t      P>|t|      [0.025      0.975]
-------------------------------------------------------------------------------------
const                 4.2239      0.009    448.912      0.000       4.205       4.242
number_of_ratings    -0.0179      0.001    -17.631      0.000      -0.020      -0.016
==============================================================================
Omnibus:                     5957.727   Durbin-Watson:                   1.919
Prob(Omnibus):                  0.000   Jarque-Bera (JB):           261509.787
Skew:                          -2.213   Prob(JB):                         0.00
Kurtosis:                      27.658   Cond. No.                         26.7
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
In [227]:
fig, axes = plt.subplots(figsize=(10,10))

plt.scatter(x, y, alpha = 0.1)

polyline = np.linspace(min(x), max(x), max(x)-min(x))
polyline2 = np.linspace(0, 5, 1000)

axes.plot(polyline, lin_model.params[0] + lin_model.params[1] * polyline,
          color="red", label="Линейная")
axes.plot(polyline, hyper_model.params[0] + hyper_model.params[1] * (1 / (polyline + 1)),
          color="yellow", label="Гипербола")
axes.plot(polyline, sqrt_model.params[0] + sqrt_model.params[1] * polyline ** 0.5,
          color="green", label="Корень")
axes.plot(polyline, log_model.params[0] + log_model.params[1] * np.log(polyline + 1),
          color="blue", label="Логарифм")
plt.xticks(fontsize=9)
plt.yticks(fontsize=9)
#plt.xlim(-10, 0.2 * 1e6)

axes.legend(loc="best", fontsize=10)
plt.show()

R:

In [228]:
%%R

x <- df2$number_of_ratings
y <- df2$average_rating

lin_model <- lm(y ~ x)
print(summary(lin_model))
hyper_model <- lm(y ~ 1/(x+1))
print(summary(hyper_model))
sqrt_model <- lm(y ~ sqrt(x))
print(summary(sqrt_model))
log_model <- lm(y ~ log(x+1))
print(summary(log_model))
Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0660 -0.1884  0.0040  0.1928  0.9340 

Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept) 4.066e+00  3.474e-03 1170.342  < 2e-16 ***
x           2.787e-08  9.799e-09    2.845  0.00445 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3352 on 9998 degrees of freedom
Multiple R-squared:  0.0008088,	Adjusted R-squared:  0.0007088 
F-statistic: 8.093 on 1 and 9998 DF,  p-value: 0.004453


Call:
lm(formula = y ~ 1/(x + 1))

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0686 -0.1886  0.0114  0.1914  0.9314 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 4.068577   0.003354    1213   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3354 on 9999 degrees of freedom


Call:
lm(formula = y ~ sqrt(x))

Residuals:
    Min      1Q  Median      3Q     Max 
-4.0755 -0.1893  0.0075  0.1911  0.9249 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.075e+00  4.267e-03 955.063  < 2e-16 ***
sqrt(x)     -3.650e-05  1.399e-05  -2.609  0.00909 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3353 on 9998 degrees of freedom
Multiple R-squared:  0.0006804,	Adjusted R-squared:  0.0005805 
F-statistic: 6.808 on 1 and 9998 DF,  p-value: 0.009091


Call:
lm(formula = y ~ log(x + 1))

Residuals:
    Min      1Q  Median      3Q     Max 
-4.2239 -0.1767  0.0214  0.1944  0.8571 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.223915   0.009409  448.91   <2e-16 ***
log(x + 1)  -0.017876   0.001014  -17.63   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.3303 on 9998 degrees of freedom
Multiple R-squared:  0.03015,	Adjusted R-squared:  0.03006 
F-statistic: 310.9 on 1 and 9998 DF,  p-value: < 2.2e-16

In [229]:
%%R -w 1000 -h 1000

plot(x, y)

abline(lin_model, col='red')

xs = seq(from = min(x), to = max(x), length.out = max(x)-min(x))
hyper_ys = predict(hyper_model, newdata = list(x=seq(from = min(x), to = max(x), length.out = max(x)-min(x))))
matlines(xs, hyper_ys, lwd = 2, col = 'yellow')

sqrt_ys = predict(sqrt_model, newdata = list(x=seq(from = min(x), to = max(x), length.out = max(x)-min(x))))
matlines(xs, sqrt_ys, lwd = 2, col = 'green')

log_ys = predict(log_model, newdata = list(x=seq(from = min(x), to = max(x), length.out = max(x)-min(x))))
matlines(xs, log_ys, lwd = 2, col = 'blue')

Полученные модели, особенно последние три, дают очень близкий результат (в частности, их остаточные стандартные ошибки крайне близки).

Также по коэффициентам регрессии видно, что количество оценок практически не влияет на рейтинг.